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An  analytical  solution  is  developed  for  studying  fluid 
flow  and  heat  transfer  in  the  end  cap  of  a  single-pass, 
return  flow  heat  exchanger,  under  stokes  conditions  where  the 
inertial  effects  are  negligible  in  comparison  with  those  o£ 
viscosity.     The  cap  has  the  configuration  of  an  oblate 
spheroid  and  the  flow  is  considered  to  be  in  a  steady  state. 
A  low  speed  flow  of  an  incompressible  medium  is  considered 
for  which  the  stream  function  is  expressed  in  a  series  of 
Legendre  polynomials  which  converge  rapidly  for  the  ranges 
of  the  parameters  studied  in  this  work. 

In  the  study  of  the  heat  transfer  in  the  end  cap,  the 
temperature  field  is  solved  using  a  parameter  perturbation 
method.     The  zeroth  order  perturbation  corresponds  to  a  pure 
conduction  problem,   and  the  convective  effect  is  accounted  for 

xii 


in  the  first  order  perturbation  solution  which  is  also  the 
highest  order  considered  in  this  work. 

The  computed  stream  functions  are  used  to  plot  stream- 
lines, which  are,   in  turn,  utilized  to  determine  the  velocity 
and  pressure  fields.     The  temperature  fields  are  used  to  plot 
isotherms  from  which  the  local  Nusselt  number  is  also  determined. 
The  analysis  reveals  that  there  is  a  point  of  singularity  for 
the  vorticity  at  the  rim  of  the  feeder  tube.     A  recirculation 
flow  is  also  found  when  the  annulus  region  has  an  area  that  is 
large  compared  to  the  feeder  tube. 
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CHAPTER  I 
INTRODUCTION 

Of  the  numerous  heat  exchangers  developed  in  the  past, 
the  coaxial,  return-flow  heat  exchanger  is  one  of  the  most 
interesting.     It  consists  of  two  coaxial  tubes,  one  inside 
the  other,   see  Figure  1.     The  round  cap  at  one  end  serves  as 
a  guide  for  directing  the  flow  from  one  passage  to  the  other. 
There  are  many  variations  in  which  this  exchanger  can  be  used. 
The  inner  tube  can  be  used  as  the  feeder  for  fluid  flow. 
Then,  the  annulus  serves  as  the  return  channel.     The  fluid 
flow  direction  may  be  reversed  using  the  annulus  as  t,he 
feeder  and  the  inner  pipe  as  the  return  passage.     In  addition, 
several  flow  patterns  may  develop  within  the  passages  depend- 
ing on  the  tube  dimensions.     For  example,   if  the  radius  oi 
the  outer  pipe  is  large  compared  with  that  of  the  inner  pipe, 
a  recirculation  flow  appears  near  the  inner  end  of  the  inner 
pipe.     This  circulation  is  not  present  for  moderate  radiiis 
ratios.     Furthermore,   if  the  inner  pipe  extends  toward  the 
end  cap  such  that  a  small  gap  is  formed  between  the  tube  and 
the  round  cap,  this  obstructs  the  flow  such  that  another  flow 
condition  prevails. 

Notice  that  those  conditions  used  to  characterize  the 
fluid  flow  and  heat  transfer  inside  tubes  and  annuli  may  also 
be  used  for  the  study  of  the  coaxial  exchanger.     Thus,  from 
the  fluid  flow  point  of  view,  there  may  be  laminar  or  turbulent 
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flow,  and  developing  and  fully  developed  flow.  On  the  other 
hand,  from  the  heat  transfer  standpoint,  there  may  be  a  pre- 
scribed temperature  or  a  heat  flux  imposed  as  boundary  condi- 
tions. The  thermal  boundary  layer  may  be  fully  developed  or 
may  be  in  the  developing  stage.  The  coaxial,  return-flow 
heat  exchanger  thus  encompasses  an  array  of  problems  for 
study . 

It  may  be  noted  that  the  coaxial  heat  exchanger  has  been 
used  by  Owens-Illinois  to  develop  the  SUNPAK  solar  collector, 
see  Figure  2.     Coaxial  heat  exchangers  are  inserted  into  a 
common  manifold  at  the  center  of  the  tube  bank  such  that  the 
fluid  flow  in  the  exchanger  is  connected  in  series.  SUNPAK 
collectors  use  an  evacuated  envelope  for  the  exchanger  to  re- 
duce the  heat  loss.     Their  thermal  efficiency  is  about  60%  [1]. 

The  coaxial  heat  exchanger  has  also  been  used  in  Alkylation 
contactors  in  which  the  exchangers  form  a  tube  bundle  erected 
in  a  vertical  position,  see  Figure  3.     One  of  the  advantages 
of  this  coaxial  exchanger  in  this  erected  installation  is 
that,  since  it  is  mounted  at  one  end,   it  is  free  to  expand  at 
the  other  end,  thus  minimizing  stress  and  seal  problems. 

While  the  exchanger  has  been  considered  useful  for  many 
applications,  there  has  been  little  work  done  to  study  its 
fluid  flow  and  heat  transfer.     This  motivated  the  research 
given  in  this  work.     This  dissertation  presents  an  analytic 
solution.     Part  of  the  motivation  for  this  approach  is  that 
numerical  solutions  (finite  difference  or  finite  element), 
while  they  can  always  be  used  for  solution,  must  be  tested 
for  convergence.     Such  a  test  may  be  costly  for  a  complicated 


4 


Figure  2.     Sunpak  Series  Collectors 
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Figure 


3. 


Field  Tubes  in  an  Alkylation 


Contactor 
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system  like  the  heat  exchanger.     It  is  necessary  to  accept 
certain  restrictions  in  order  to  achieve  an  exact  solution. 
However,  it  does  not  suffer  from  the  convergence.  Notice 
also  that  the  fluid  flow  and  heat  transfer  inside  tubes  and 
annuli  have  been  thoroughly  studied  and  well  documented  in 
the  literature,  only  the  end  cap  poses  as  a  new  problem  for 
analysis.     The  exact  solution  is  thus  dedicated  to  that,  and 
this  also  defines  the  scope  of  this  work. 

In  the  chapters  that  follow,   a  literature  survey  will  be 
provided  in  Chapter  II.     The  formulation  of  the  problem  and 
the  procedure  for  solution  that  leads  to  the  determination 
of  the  stream  function,   and  consequently  the  velocity  and 
the  pressure  distributions  in  the  cap,   is  presented  in 
Chapter  III.     In  Chapter  IV  an  analysis  of  the  heat  transfer 
problem  is  provided  for  the  solution  of  low  Prandtl  n\amber 
situations.     The  results  and  discussions  are  summarized  in 
Chapter  V.     Finally,  some  conclusions  and  recommendations  are 
made  in  Chapter  VI. 


CHAPTER  II 
LITERATURE  REVIEW 


The  coaxial,   return-flow  heat  exchanger  may  be  struc- 
turally divided  into  three  components:     the  inner  feed  pipe, 
the  annulus,  and  the  end  cap.     Studies  of  fluid  flow  and  heat 
transfer  in  these  three  components  will  be  briefly  reviewed 
in  this  chapter. 

It  is  noted  that  the  inner  pipe  and  its  surrounding 
annulus  are  two  components  in  which  the  fluid  flow  and  heat 
transfer  problems  have  been  thoroughly  studied  in  the  litera- 
ture.    Review  monographs  are  also  available  for  these  systems. 
Hence,  their  review  given  here  will  be  brief,  and  their  inclu- 
sion is  mainly  for  the  sake  of  completeness.     For  an  orderly 
presentation  of  the  materials  that  follow,  those  monographs 
will  be  cited  first.     This  will  be  followed  by  a  brief  review 
of  the  important  papers  published  dealing  with  these  subjects. 

Fluid  Flow  and  Heat  Transfer  in  Tubes 
Fluid  flow  and  heat  transfer  in  tubes  and  ducts  are  a 
topic  of  typical  textbooks,   e.g.    [2-5] .     Review  of  studies 
of  fluid  flow  and  heat  transfer  in  ducts  may  be  found  in  a 
series  of  monographs  [6] .     The  most  comprehensive  review  was 
published  by  Shah  and  London  [5j  in  1978. 

A  review  of  the  literature  reveals  that  the  first  analyti- 
cal solution  for  heat  transfer  in  a  closed  conduit  was  prescribed 
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by  Graetz  [7]  in  1883.     Nusselt  made  an  independent  analytic 
solution  in  1910   [8].     Dryden  et  al.    [9]  compiled  solutions 
for  fully  developed  laminar  flow  inside  ducts  of  various 
geometries,   in  1932. 

Three  methods  of  analysis  have  been  applied  to  the  flow 
in  the  entrance  region.     One  method  applies  the  integral  form 
of  the  equations  of  momentum  and  mass  continuity  to  the 
boundary  layer  as  it  develops  along  the  tube  wall.     The  velo- 
city profile  is  approximated  by  a  suitable  polynomial  as  in  the 
Karman-Pohlhausen  method. 

The  second  method  subdivides  the  entrance  region  into 
two  zones.     In  the  zone  near  the  tube  entrance,   a  boundary- 
layer  equation  is  used.     An  approximate  solution  is  obtained 
by  means  of  a  perturbation  series  which  takes  the  fully  devel- 
oped velocity  distribution  as  a  limiting  case  in  the  second 
region  where  the  flow  moves  far  downstream  from  the  entrance 
section.     The  flow  development  throughout  the  entrance  region 
is  found  by  patching  together  the  upstream  and  downstream 
solutions  at  some  intermediate  location. 

In  the  third  method,   a  solution  for  the  velocity  is  ob- 
tained that  is  continuous  radially  and  axially  from  the  en- 
trance to  the  fully  developed  region.     Langhaar  [10]  in  1942 
linearized  the  Navier-Stokes  equations  and  solved  to  obtain  a 
family  of  velocity  profiles  in  the  form  of  Bessel  functions. 
The  parameter  of  this  family  was  tabulated  against  the  axial 
coordinate  in  dimensionless  form,   and  ac  a  result,  the  transi- 
tion length  was  also  obtained.     Sparrow  et  al .    [11]  linearized 
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the  momentum  equation  using  the  mean  velocity  along  the  flow 
direction.     Sparrow's  solution  appears  to  be  more  general 
since  it  is  applicable  to  both  a  circular  pipe  and  a  parallel- 
plate  channel.     A  further  improvement  was  made  by  Fleming 
and  Sparrow  [12]  in  1969  allowing  application  to  the  entrance 
region  of  ducts  of  arbitrary  cross  section. 

Analysis  for  heat  transfer  in  a  circular  tube  with  fully 
developed  velocity  and  temperature  profiles  is  straight  for- 
ward when  the  wall  temperature  or  the  heat  flux  is  uniform. 
For  an  arbitrary  circumferential  heat  flux  but  constant  axial 
heat  input,  a  superposition  method  can  be  used  for  solution 
[13].     When  the  flow  is  fully  developed  but  heat  transfer  is 
in  the  developing  stage,  the  problem  can  be  solved  by  separa- 
tion of  variables.     A  specific  example  of  this  type  is  the 
case  of  uniform  surface  temperature  when  axial  heat  condution 
can  be  neglected;  the  solution  in  this  case  is  in  the  form  of 
Graetz  functions.     The  eigenvalues  for  the  infinite  series 
in  this  solution  can  be  obtained  by  numerical  approxima- 
tion.    These  methods  have  also  been  attempted  by  Sellars 
et  al.   [14]  in  the  study  of  similar  problems  with  uniform 
heat  flux  as  boundary  conditions. 

The  thermal  entry  length  solutions  described  above  form 
building  blocks  for  studying  heat  transfer  in  pipes  with  sur- 
face temperature  variations.     Duhumel's  theorem  is  useful 
in  this  instance.     A  similar  analysis  can  be  used  to  study 
problems  with  axial  heat  flux  variations.     However,  the 
problem  becomes  more  complicated  when  both  velocity  and 
temperature  profiles  develop  simultaneously.     The  Prandtl 
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number  plays  a  leading  role  in  this  problem,   and  numerical 
methods  are  often  used  for  solution. 

Fluid  Flow  and  Heat  Transfer  in  Annuli 
Murakawa  [15]  solved  the  continuity  and  momentum  equa- 
tions for  the  velocity  distribution,  pressure  drop,   and  entry 
length  for  developing  flow  inside  a  circular  annulus.  An 
approximate  solution  was  also  developed  by  him  for  studying 
heat  transfer  in  the  thermal  entry  length.     Hatton  and  Quarmby 
[16]  developed  a  solution  to  the  Graetz  problem  in  an  annulus 
with  a  heated  core  and  an  insulated  outer  wall.     Both  uniform 
temperature  and  uniform  heat  input  on  the  inner  wall  were  con- 
sidered for  a  wide  range  of  arbitrarily  chosen  radius  ratios. 
Their  solution  was  also  extended  to  the  case  of  axial  variation 
of  the  core  heat  input  and  to  the  case  in  which  heat  was  trans- 
ferred at  uniform  but  unequal  rates  at  both  walls. 

The  most  comprehensive  study  of  heat  transfer  in  an  an- 
nulus was  made  by  the  group  led  by  Kays  and  Reynolds  at 
Standford.     A  series  of  four  papers  was  published  by  them,  with 
the  first  one  including  a  bibliography  covering  those  annuli 
papaers  published  before  1961.     Four  classes  of  problems  were 
formulated  in  this  first  paper  [17],  which  consisted  of  (i)  a 
step  temperature  on  one  wall,  with  the  other  wall  maintained 
at  the  inlet  temperature,   Cii)  ^  step  heat  flux  from  one  wall, 
with  the  other  wall  kept  at  an  adiabatic  condition,   (iii)  a 
step  temperature  on  one  wall,  with  the  other  wall  insulated, 
and  finally  (iv)  a  step  heat  flux  from  one  wall,  with  the 
other  maintained  at  the  inlet  temperature.     The  utility  of 
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these  four  basic  problems  was  also  described  in  this  paper 
along  with  a  superposition  technique  useful  for  the  solution 
of  more  complicated  problems.     This  paper  was  concluded  with 
a  brief  description  of  the  experiments  which  they  used  to 
verify  analysis. 

In  the  second  paper  [18],  solutions  for  those  four  basic 
problems  were  developed  for  the  case  of  a  fully  developed  flow. 
In  the  analysis,  the  viscous  generation  of  heat  and  the  axial 
heat  conduction  effect  were  neglected  in  the  energy  equation. 
The  temperature  distributions,  obtained  from  the  solution  of 
eigenvalue  problems,  were  verified  with  experiments.  Super- 
position of  solutions  of  basic  problems  was  also  illustrated 
by  means  of  examples. 

In  the  third  paper  [19],  heat  transfer  with  turbulent 
flow  was  considered.     The  flow  was  taken  .tr>  be  fully  developed 
inside  an  axi symmetrically  heated  annulus.     Fluid  properties 
were  treated  as  constant,   and  axial  heat  conduction  was  again 
neglected  in  the  analysis.     Asymptotic  solutions  (fully  developed 
velocity  and  temperature  profiles)  were  developed  for  a  wide 
range  of  radius  ratios,  Reynolds  numbers,   and  Prandtl  numbers. 
The  analysis  was  made  on  the  basis  of  empirical  velocity  and 
eddy  diffusivity  profiles,   and  experiments  were  carried  out 
to  verify  the  analysis  at  Pr=0.7.     A  superposition  technique 
was  also  demonstrated  for  solving  problems  with  axisymmetric 
heating  from  both  walls  of  the  annulus. 

Heat  transfer  for  laminar  flow  in  an  annulus  with 
simultaneously  developing  velocity  and  temperature  profiles 
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was  studied  by  Heaton  et  al.    [20]  in  the  final  paper  in  the 
series.     A  solution  was  obtained  first  for  the  hydrodynamic 
problem,  and  then  for  the  combined  hydrodynamic  and  thermal 
problem  with  an  integral  method.     Excellent  agreement  was 
again  found  between  the  experiment  and  the  analysis. 

As  a  further  extension  to  what  was  published  earlier 
in  the  series,  Sutherland  and  Kays  [2l]  studies  heat  transfer 
in  an  annulus  with  variable  circiimf erential  heat  flux  but 
constant  heat  rate  along  the  tube  axis.     The  velocity  and 
temperature  profiles  were  fully  developed,  and  the  analysis 
was  carried  out  for  both  laminar  and  turbulent  flows. 

It  is  noted  that,   for  those  studies  cited  above,  the 
fluid  properties  have  been  treated  as  constant.  Shumway 
and  McEligot  [22]  studied  heat  transfer  in  an  annulus  with 
temperature  dependnent  transport  properties.     Numerical  solu- 
tions were  carried  out  for  both  fully  developed  and  uniform 
entry  velocity  profiles. 

Fluid  Flow  and  Heat  Transfer  in  End  Caps 
Unlike  the  pipe  and  the  annulus,  there  has  been  little 
attention  given  to  the  study  of  fluid  flow  and  heat  transfer 
in  the  end  cap.     There  is  a  lack  of  literature  published 
dealing  directly  with  the  fluid  flow  problem.     However,  a  re- 
lated problem  of  magnetohydrodynamic  flow  cap  has  been  in- 
vestigated [23,24].     In  this  magnetohydrodynamic  problem, 
the  fluid  is  set  in  motion  due  to  the  Lorentz  force  resulting 
from  the  electric  current  and  the  associated  magnetic  field. 
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The' case  of  flow  field  induced  in  an  incompressible  viscous 
fluid  was  studied  for  discharging  either  from  a  point  source 
[25]  or  a  plane  electrode  [26]  located  at  the  center  of  the 
plane  side  of  the  hemisphere.     For  the  point  source  problem, 
an  analytic  solution  has  been  developed  for  the  slow  viscous 
flow  and  a  numerical  solution  for  the  nonlinear  problem.  As 
one  might  expect,  because  of  the  axisymmetry  of  the  position 
of  the  electrode,  stream  lines  in  an  axial  cross  section  form 
two  closed  loops.     For  the  plane  electrode  case,  the  vorticity 
has  a  singularity  at  the  rim  of  the  electrode.     When  the  width 
of  the  annulus,  that  is  the  remaining  of  the  plane  side  of 
the  hemisphere,  is  small  relative  to  that  of  the  electrode, 
an  eddy  is  formed  about  the  rim  of  the  electrode.     This  eddy 
diminishes  in  size  as  the  annulus  width  is  increased. 

Fluid  Flow  and  Heat  Transfer  in  the  Total  System 
There  have  been  few  studies  published  for  fluid  flow  and 
heat  transfer  in  the  total  system  of  the  coaxial  exchanger. 
Hurd  [27]  studied  the  mean  temperature  difference  for  the 
exchanger  used  in  the  Alkylation  contactor  in  which  the  inner 
feed  pipe,  or  the  annulus,  was  used  to  deliver  the  coolant 
and  was  thus  at  the  shell  side  of  the  contactor.     Based  on 
the  relative  directions  of  flows  between  the  annulus  fluid 
and  its  shell-side,  heat  exchanging  fluid,  four  flow  conditions 
may  arise.     Studies  were  made  by  Hurd  to  develop  the  mean  temp- 
erature difference  equations  and  curves  for  thse  four  condi- 
tionss.     Hurd's  paper  centers  on  the  study  of  the  application 
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aspect  of  the  heat  exchanger.  It  provides  little  insight  into 
the  mechanisms  of  fluid  flow  and  heat  transfer. 

More  recently  Moan  [28]  and  Beekley  and  Mather  [l]  pub- 
lished performance  data  of  coaxial  exchangers  fitted  with 
evacuated  envelopes  to  reduce  the  heat  loss.     The  exchanger 
was  investigated  for  solar  energy  applications,   and  the  analy- 
sis presented  was  focused  on  the  total  system  behavior. 


CHAPTER  III 
ANALYSIS  OF  THE  FLUID  FLOW  PROBLEM 


Governing  Equations 


Consider  a  steady,  incompressible  flow  of  a  viscous 
fluid  ejecting  from  a  circular  tube  of  radius  d  into  an 
oblate  hemispheroidal  cap  and  delivered  through  an  annular 
passage  of  outer  radius  a,  see  Figure  4.     The  equitorial  and 
polar  radii  of  the  spheroid  are  a  and  c,   respectively.  The 
axes  of  tubes  lie  on  the  axis  of  the  cap,  which  is  also  the 
axis  of  symmetry  of  the  heat  exchanger. 

Let  V  be  the  fluid  velocity.     Under  the  condition  of 
constant  properties,  no  body  force,   and  sufficiently  slow 
motion  of  the  fluid  such  that  the  squares  and  products  of  the 
velocity  components  are  negligible,   the  continuity  equation 
and  the  Navier-Stokes  equations  can  be  simplified  to  be 


where  p,   p,   and  v  have  their  usual  meaning  of  pressure, 
density,   and  kinematic  viscosity,  respectively. 

It  is  convenient  to  make  use  of  the  oblate  spheroidal 
coordinates  (C,  u,  (}>)  which  are  related  to  the  cylindrical 
coordinates  (r,  z,   (jj)  by 


V 


V  = 


0 


(3.1) 


-Vp  +  pvV^V  =  0 


(3.2) 


r  =  d(l  +  c2)^(l-y2)4 


(3.3) 
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z  =  dcy  (3.4) 

Here  the  origin  oi  the  coordinate  system  (r,  z,  ())  )  is  taken 
to  be  the  center  of  the  spheroid,  and  the  positive  z-axis  is 
pointed  toward  the  cavity  side  of  the  cap  and  is  oriented 
along  the  axis  of  symmetry.     Equations  3.3  and  3.4  are  derived 
in  Appendix  A. 

The  inlet  region  can  now  be  expressed  as 

T  <_d  and  z  =  0 

in  the  cylindrical  system,  or 

Oj<  y   £  1  and  C  =  0 

in  the  spheroidal  system.     Similarly,  the  annular  passage  is 
in  the  region 

d  _<  r  j<  a  and  z  =  0 

or 


y  =  0  and  0  <  ^  < 

—  ^  —  ^o 


The  surface  of  the  cap  is  given  byC  =  ^Q>0£l^i.^' 

fluid  region  byOj<y£l,   0  <  i.  <  r,^,   see  Figure  4.     With  the 

geometry  given,  the  equatorial  and  polar  radii  of  the  spheroid 

can  be  related  to  d  and  c  by 

o 

a  =  d(l+  C^)*  and  c  =  d? 
^       "^o  o 
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Figure  4.     Section  View  of  the  Oblate  Spheroid 
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Since  the  velocity  field  V  is  axisymmetric  and  thus  two 
dimensional,   it  can  be  expressed  in  terms  of  the  stream  func- 
tion which  satisfies  the  continuity  equation  identically, 
The  velocity  components  u  and  v  can  be  derived  (see  Appendix  A) 
as 


V  = 


d"(?'  +  li") 


2  N  t 


-1 


2^* 


(1+CO 


9y 


(i-y^)* 


_3Y 

9? 


(3.5) 


Taking  curl  of  Equation  3.2,   and  making  use  of  vector  identities 
V^V=V( V'V)-Vx(VxV)  and  VxVp=0  and  also  the  continuity  equation 
3.1,   it  can  be  derived  that 


or 


VxVx(VxV)  =  0 


VxVxw  =  0 


(3.6) 


where  u)  is  vorticity,  defined  as 


w  =  VxV 


(3.7) 


Finally  introducing  ^=vd<l')  and  C=?qX  and  making  use  of  Equation 
3.5,  one  can  transform  Equation  3.6  to 


(x^+y^C;2) 


D^^'=  0 


(3.8) 


where 


(x2+;-2) 


9'     .  2^  9^ 


3x' 


+  (i-y^) 


9y' 
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Equation  3.8  becomes  the  transformed  governing  equation 
for  the  fluid  flow  problem. 

Boundary  Conditions 

The  governing  Equation  3.8  will  be  solved  subjected  to 
a  set  of  boundary  conditions  described  below. 

Along  the  axis  of  symmetry  (i.e.   y=l),  the  stream  func- 
tion must  have  a  constant  value.     It  is  convenient  to  take 
this  value  to  be  zero.     Physically,   this  characterizes  the 
condition  that  the  velocity  component  v,   cross-flow  is  zero 
along  this  axis.     At  the  spheroidal  surface  (x=l),  refer  to 
Figure  5,  the  velocity  components  u  and  v  both  vanish.  Be- 
sides, this  surface  can  be  considered  as  an  extension  of  the 
axis  of  symmetry  (y=l)  from  the  stream  function  standpoint. 
Hence,  ¥  is  zero  on  this  surface. 

For  the  low  Reynolds  number  flow  that  is  of  interest  in  . 
this  work,  the  fluid  velocity  emerging  from  the  feed  pipe  and 
exiting  to  the  annulus  region  is  expected  to  be  nearly  parallel 
to  the  axis  of  symmetry.     It  is  common  practice  [29,30]  to 
consider  the  velocity  component  along  both  x=0  and  y=0  planes 
to  be  negligibly  small,  refer  to  Figure  5.     Finally,   the  fluid 
velocity  emerging  from  the  feed  pipe  is  assumed  to  be  known 
[31,32].     The  stream  function  can  thus  be  determined  at  x=0. 

These  conditions  can  be  summarized  as  follows 


¥(x,y=l)  =  0 
1'Cx=l,u)  =  0 


(axis  of  symmetry) 
(spheroidal  surface) 


(3.9) 
C3.10) 
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Figure  5.     Boundary  Conditions  for  the  Fluid  Flow  Problem 
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m^i^  =  o  (3.11) 

oX 


^'^^9x^'^''  "  Q  (inlet)  (3.12) 


H'(x=0,y)  =  f(y)  (3.13) 


=  0  (exit)  (3.14) 


They  can  also  be  put  into  perspective  in  Figure  5. 

Development  of  a  Solution  Technique 
The  problem  at  hand  can  be  solved  by  expanding  the 
stream  function  in  a  series  of  Legendre  polynomials 


^(x.y)  =  (l-y^       E     g2n_i(^)  ^2n-l^^^  (3.15) 

n=l 


Each  term  in  this  series  satisfies  Equations  3.9  and  3.14 
identically.     To  satisfy  Equations  3.10  through  3.12  calls  for 
meeting  new  conditions  as  follows 

g2n_iCl)  =  0  (3.16) 
g^„_l(l)  =  0  (3.17) 
S2n-1^°^  =  0  (3.18) 

As  mentioned  earlier,  Equation  3.13  signifies  that  the 
velocity  is  known  at  the  exit  of  the  feed  pipe.     It  is  possible 
to  expand  f(y)  in  a  series  of  ( 1-y  ^  )P2j^_;[^(  y )  over  the  region 
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0  <  y  <  1  provided  that  the  function  satisfies  Dirichlet 
conditions  in  the  interval  -1  1  y  1  1  [33].     For  the  coaxial 
heat  exchanger  problem  under  investigation,  let  the  velocity 
emerging  from  the  feed  pipe  be  fully  developed.     Then,  the 
function  f(y)  can  be  related  to  y  (see  Appendix  B)  as 

¥(0,y)  =  f(y)  =  iReCl-y")  (3.19) 

where  the  Reynolds  number  is  defined  as  2Vd/v,  and  V  refers 
to  the  mean  velocity  at  the  pipe  exit . 

In  Equation  3.19,  f(y)  can  be  expanded  in  a  series 

following  the  pattern  given  in  Equation  3.15.  Thus, 

'l'(0,y)  =  (iRe)(l-yM   [(3/5  P^Cy)  +  1/15  p'^iM)]  (3.20) 

A  comparison  of  this  equation  with  Equation  3.15  at  x=0 
provides  the  following  conditions 

g-^(O)  =  3/10  Re  (3.21) 

g3(0)  =  1/30  Re  (.3.22) 

and 

g2^_-^(0)  =0  for  n  ^  3  C3.23) 

These  will  be  used  as  boundary  conditions  together  with 
Equation  3.16  through  Equation  3.18  in  the  solution  for  g. 
Notice  that  the  function  itself  will  satisfy  a  set  of  ordinary 
differential  equations  whose  derivation  follows  next. 

The  stream  function  ^  should  satisfy  Equation  3.8.  Hence, 
the  series  of  the  Legendre  polynomials  3.15  is  substituted  into 
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Equation  3.8,  and  after  some  manipulation  (see  Appendix  C), 
giving 


n=l 


■*-CG2n-l(")^'P2n-l(^)>  "  ^4^" ( 1-y ^  )P2^_,(  y )  ] ' 


+  2[(n+l)(2n-3)x2-2c-^]P2^_^(y) 


+2[(n+l)(2n-3)C-2+6?-2]y2p^^_^(^)}G2^_^(x)  =  0  (3.24) 


where 


Sn-l^'')  =        ""^o  )^2n-l^'')"2''^2''~^^^2n-l^''^  (3.25) 

In  order  to  determine  G  from  this  equation,   it  is  appro- 

priate  to  expand  terms  involving  V^^^2n-1^^^  y[(         )P2n_i^^^-^  ' 

in  terms  of  Legendre  polynomials.     Use  is  thus  made  of  two 
identities  [33] : 

u2p'       (   )  ^    2n(2n-l)   '  (2n-l)  +  (2n-2)(4n+l)  ■ 

^  ^2n-l^'^''       (4n-l)(4n+l)  ^2n+l*^^''  (4n+l)(4n-3)  ^2n-l'-^^ 

^  (4MffiE^P;n-3(u)  (3.26) 

y[(l-y')P2^_^(y)]  '  = 

2n(2n-l)  r2n-l     '                        4n-l  .    .       2n     „ '       ,    , -, 

4n-l       '-4n-3  '2n-3^^^  "  (4n+l)(4n-3)  ^2n-l'^^ 4H+I  ^2n+l^^-*-' 

(3.27) 
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which  are  substituted  into  Equation  3.24,     The  series  in  the 
resulting  equation  is  then  expanded,   and  finally  those  terms 
containing  derivatives  of  the  Legendre  polynomials  of  the 
same  degree  (i.e.  terms  of  P^,  P^,  etc.    .    .    .)  are  grouped 
together  and  set  to  zero.     A  general  formula  is  derived, 


^(2)  Qu        +  JO)              ^  (2)              ^     (1)  • 

^2n-3^2n-3  ^2n-3^2n-3  °2n-1^2n-l  °2n-1^2n-l 

+  b^°\G_     ,  +  c^^],G"  +  c^°2iGo  o-i  =  0  (3.28) 

2n-l  2n-l        2n+l  2n+l  2n+l  2n+l 


In  this  equation,   coefficients  a,  b  and  c  have  been  identified 

with  both  superscripts  and  subscripts  to  show  where  they  belong. 

The  superscripts  are  specifically  used  to  relate  the  order  of 

the  derivatives  of  the  G  function  with  which  they  are  multiplied. 

These  coefficients  are  functions  of  ^  , n  and  x, 

o'  ' 


_C2)     ^  fx^  +  C-M  2(n-l)(2n-3)       2  , 
^2n-3  ^o  ^   C4n-3)(4n-5)   ^o  {d.Zii) 


(0)     _  4(n-l)(2n-3)(-2n^+9n-10)       2  .0  nr.. 

^2n-3  (4n-3)(4n-5)  L-^.^V) 


b(2)     =  (x^  +  L-^)(x^+  (8n^-4n-3)   2x  ,0  on n 

^2n-l  ^o  C4n+l)(4n-3)  ^o  ^  (d.Jl) 


^2^-1  "  -4x(x2+c-2)  (3.32) 


^(0)     ^  8n(2n-l)-2(2n^-n+3)(8n^-4n-3)  2 
2n-l  C4n+l)(4n-3)  ^o 


-2[(n+l)(2n-3)x2-2c-^]  (3.33) 
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p(2)     -  rx=^  +  r-M  2(2n+l)(n+l)       2  .3 
*'2n+l  "  (4n+3)(4n+l)  ^-  •''^^ 


(0)     _  -4(n+l)(n+2)(2n  +  l)(2n+3)       2  .3 
''2n+l  "  (4n+3)(4n+l)  {d.cib) 


and  a^^^  and  a^^^  are  taken  to  be  zero. 

Equation  3.28  forms  an  infinite  set  of  second  order, 
linear  differential  equations;  but,   it  was  g  conditions  that 
were  derived  earlier.     Hence,  Equation  3.28  is  changed  to  a 
differential  equation  in  g  by  using  the  definition  of  G, 
Equation  3.25.     This  leads  to  a  fourth  order,   linear  differen- 
tial equation  of  g  as 


^2n-1^2n-l      "2n-1^2n-l      ^2n-1^2n-l  ^2n-1^2n-l 

^  ^2n-1^2n-l  =  -^2n-l  "  ^2n-l  ^"=1'  2.    .    .    .    ,  N)  (3.36) 


where 


=  0  (3.37) 

2n-l        2n-1^2n-3        2n-1^2n-3  '^2n-1^2n-3 

^  Sn-lS2n-3  3,    .    .    .    ,  N)  (3.38) 


^2n-l       ^2n-1^2n+l      ™2n-1^2n+l  ^2n-1^2n+l 

^  ^2n-1^2n+l  2,    .    .    .    ,  N-1)  (3.39) 


S2N-I  =  0 


(3.40) 
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In  the  above  formulation,  the  series  is  assumed  to  terminate 

at  N,     Also  notice  that  the  new  coefficients  A,  B,    .    .    ,  are 

related  to  the  a,  b,   c  coefficients  mentioned  earlier  (Equations 

3,29  through  3,35).     The  detailed  derivation  for  these  equations 

is  given  in  Appendix  C, 

Several  features  in  Equations  3.36  through  3,40  are 

worthy  of  note.  t   and  S_     ,   appear  on  the  right  hand 

^  2n-l  2n-l 

side  of  Equation  3.36.     Since  ^211-1  ''"^  related  to  g2n-3 

through  Equation  3.38  and  S2^_-|^  to  S2^+i  through  Equation  3.39 

and  the  fact  that  the  subscript  n  here  designates  the  nth 

term  in  the  series  expansion,  these  g  terms  for  succeeding 

expansions  in  the  series  are  actually  related  according  to 

Equation  3.36.     Second,   the  derivations  of  g  can  be  easily 

changed  to  finite  differences  in  the  numerical  solution. 

Equations  3.36  through  3.40  can  then  be  expressed  in 

terms  of  g(x) ,   and  these  g(x)  terms  are,  in  turn,  related  to 

their  neighboring  g  through  these  equations.     It  is  now  clear 

that  these  finite  difference  equations  must  be  solved  for  the 

g  function  for  each  n  in  the  series  expansion,  while  these  g 

functions  for  nth  expansion  are  again  related  to  their  lower 

and  higher  order  expanded  terms.     In  this  regard,  one  must  be 

aware  of  the  fact  that,   since  g„     ,   is  related  to  g^     o  and 

^2n-l  ^2n-3 

^2n+l'  ^  tri-diagonal  type  matrix  will  result.     On  the  other 
hand,  since  the  highest  order  differentials  in  Equations  3.36, 
3.38  and  3.39  are  fourth  order,   g  for  each  n  is  governed  by  a 
penta-diagonal  matrix.     With  this  in  mind,   a  scheme  for  solu- 
tion can  be  developed  as  follows  (the  details  are  again  rele- 
gated to  an  Appendix  D) . 
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As  the  first  step,  the  derivatives  of  g  are  expressed  in 
finite  difference  forms,  which  are,  in  turn,  substituted  into 
Equations  3.36,  3.38  and  3.39.  Using  h  to  represent  the  step 
size  in  the  x  direction  and  j  subscript  to  designate  position 
in  the  x  direction  change  these  equations  to  the  following 


(V-l,J-^2n-l.j  l>S2n-l,j-2       ^ -^A2n-1 ,  j^^2n-l , 

3 

^  Sn-l,j  ^'-°2n-l,j  I  )S2n-l.J-l  ^  (6A2n-l ,  j'^Sn-l ,  J 

^  ^2n-l,J  ^^)^2n-l,3  ^  ^-^A2n-1 ,  j-B2n-l ,  J  ^^Sn-l,J 
^3  ^ 
^2n-l,j  2  ^^2n-l,j  +  l       ^  ^2n-l ,  j'^^2n-l ,  j  2^^2n-l,j+2 

=  -^2n-l,J         -  S2n-l,J  (3.41) 


(^2n-l,  j"-^2n-l,  j  2^^2n-3,j-2       ("^^2n-l,j  ^■^2n-l,j  ^ 

'  J2n-l,j  ^')^2n-3,j-l  ■^(^«2n-l ,  J'^  J2n-1 ,  j         ^K2n-l,j  h^)S2n-3,J 

^-^«2n-l,j  =  ^2n-l,J  ^  ^'^2n-l,j  h^)S2n-3J  +  l 

^  («2n-l,/^2n-l,d  |)^2n-3,j-.2  =  ^2n-l,J  (3.42) 


(Wl,j-^Wl,j  |)S2n+l,j-2  ■^(-^4n-l.j^^2n-l,j  ^ 
^  (-4L2n-l,j-^Wl,j  ^-^N2n-l,J  h')S2n.l,j.l 

^  (4n-l,J^^2n-l,j  |)g2n.l,J.2  =  S2n-l,J  (^.43) 
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These  equations  can  be  slightly  simplified  by  using  boundary 
conditions  3.17  and  3.18.     Consequently,  for  Equation  3.41, 
a  penta-diagonal  matrix  results, 


'^2n-l>     =  ( 


4 

-^2n-l,2        -^2n-l,2  ^ 


4  ': 

~^2n-l,M-l  ^  ~^2n-l,M-l  ^ 


(3.44) 


where  n=l,   2,    ,    .    . ,N,   and  [t]  is  used  to  represent  the  penta- 


diagonal  matrix.     {g„     -,  }  denotes  a  vector  whose  components 

^  n "~  X 


^2n-l,l  '  ^2n-l  ' 


,gr 


2n-l,M  ' 

refers  to  the  index  for  the  end  point,  x=l 


The  subscript  M  here 


In  a  similar  fashion.  Equation  3.42  can  be  expressed  as 


^2n-l,2  ^ 


^2n-l,M-l  ^ 


0 


(3.45) 


where  n=2 ,  3, 


,  N,   and  for  Equation  3.43, 


^2n-l,2  ^ 
^2n-l,M-l  ^ 


0 


(3.46) 
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where  n=l,  2,   .    .    . ,N-1. 

These  three  matrix  equations  can  be  combined  to  give 


0 


where 


and 


0  [r3„-q 

—  n=o 


0 


—  n-o 


r  "N 

g 


'1,1 


g 


1,2 


^1,M 


g 


'3,1 


g 


3,2 


0  . 


—  n-o 


0 


0 


(3.48) 


(3.49) 


In  those  matrices  given  above,    [T] ,   [R] ,   and  [S]   are  MxM 
matrices,   and  they  are  used  in  Equation  3.47  to  construct 
the  coefficient  matrix  which  is  itself  an  NxN  matrix. 

It  is  noted  that  the  remaining  boundary  conditions  to 
be  satisfied  are  Equations  3.16,   3.21,   3.22  and  3.23.  These 
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equations  are  rewritten  by  introducing  subscripts  for  g  and 

summarized  together  with  Equations  3.17  and  3.18  in  Table  1. 

A  close  examination  of  this  table  reveals  something  interesting. 

As  always  in  this  work,  n  (in  the  first  column)  designates 
the  nth  term  in  the  series  expansion.     In  the  table,  boundary 
conditions  are  given  at  both  x=0  and  x=l,  that  is,   index  1  and 
M,  respectively.     Four  conditions  are  given  for  each  n,  which 
are  sufficient  to  solve  the  fourth  order  ordinary  differential 
equation  for  g  (Equation  3.36).     However,   as  mentioned  before, 
because  of  the  appearance  of  R  and  S  terms  on  the  right  hand 
side  of  this  equation,   g^  ^  is  also  related  to  g^  ^.     By  the 
same  token,   g„  ,  is  related  to  g,   -.   and  g^  ,,   and  g^  ,  to  g_  ^ 
and  g^  ^.     Now,  since  g  and  g'   for  n  >^  3  all  vanish  at  x=0 
and  x=l,  it  is  expected  that  gy(x)  will  be  zero  everywhere. 
Furthermore,  since  g^  -j^  is  one  ninth  of  that  g-j^      and  g^ 
is  zero,  g2(x)  is  expected  to  be  smaller  than  g-j^(x)  and  Sc^(^) 
will  be  smaller  still.     That  is  to  say,  the  Legendre  polynomial 
expansion,  Equation  3.15,  will  converge  very  rapidly  with  only 
three  terms,   at  most,  significant  in  the  series  expansion. 
Convergence  of  the  solution  is  thus  assured. 

Special  Treatment  Near  the  Rim  of  the  Feeder  Tube 
For  flow  near  the  rim  of  the  feeder  tube,   it  is  convenient 
to  consider  another  coordinate  system  (R,a)  shown  in  Figure  6. 
Writing  equations  for  c  and  n  as  functions  of  R  and  a  and 
subsequently  expanding  the  resulting  equations  using  Taylor's 
expansions  permit  writing  (see  Appendix  E) 
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Table 

1.     A  summary  of  boundary 

conditions  for  g(x) 

n 

2n-l 

conditions  given  at  x=0 

conditions  given  at  x=l 

1 

1 

g^^^=3/10  Re,  1=0 

2 

3 

g3^^=l/30  Re,   gg  1=0 

3 

5 

^5,1=0'  ^5,1=° 

4 

7 

^7,m"°'  ^7,m"*^ 

N 

2N-1 

^2N-l,l"°'  ^2N-l,l"*^ 

^2N-1,m"°'  ^2N-1,m"° 
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Figure  6.     Coordinate  System  (R,a)  at  the  Rim  of 
the  Feeder  Tube 
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C=(|)*[l+l/4(|)+l/32(|)^(l-2sina)+. . . ]  (sina/2+cosa/2)  (3.50) 
y=(f )*[l-l/4(|)+l/32(|)^(l+2sina)+. . . ]  ( cosa/2-sina/2 )  (3.51) 


for  small  R, 

R  -  idU^  +  ]i^)  (3.52) 

It  is  noted  that  the  vorticity  w  is  related  to  i>  (see 

2 

Equation  A.  7).     Since  ^=-^/vd,   and  D     is  defined  as  Equation 
A.  11.      (jj    can  be  expressed  as 


2 

vd  D  1*   ^  cQ\ 

CO  =  —5  T  1   9a  (3.53) 

d3(l+?2)*(l-li2)-(c2+yM  * 


for  small  ^  and  y, 


vd  D^H'  : 

0)  -  —5   9  (3.54) 

d^c'+y') 


Upon  using  Equation  3.52, 


(3.55) 

d  2R 


This  equation  will  be  used  later  to  show  the  singularity  for 
w  at  R=0. 

It  can  be  proved  that,   for  small  c,  (or  x)  and  y ,  the 
following  expressions  and  their  linear  combinations  are  solu- 
tions to  Equation  3.8: 
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l,y,  x,yx,CCQx2-y2),}ixSu^x,y\x^u2xS(Cjx^y-xy^),  .  .  . 

which  are  arranged  in  decreasing  order  of  magnitude.  If 
fourth  order  term  is  neglected  in  this  array,   it  can  be  shown 
that  one  that  is  most  dominant  while  satisfying  Equations 
3.12,  3.14,   ^(0 , y )=constant  (coming  from  Equation  3.19  when 
y  is  small),  and  *i'(x, 0)^1^0  is  x^.     Thus,  near  the  rim, 

4'  =  Ax^+B  (3.56) 

where  A  and  B  are  constants,  which  cannot  be  determined  with- 
out reference  to  the  overall  solution  for  the  stream  function. 
Notice  that  this  ^  is  compatible  with  the  assumed  series 
solution  for  ¥  (Equation  3.15).     As  will  be  shown  later,  g-j^(x) 
has  a  point  of  inflection,  thus  is  of  the  type  x\  Equation 
3.15  can  be  reduced  to  Equation  3.56  for  n=l  under  the  condi- 
tion of  small  y. 

Near  the  rim  of  the  feeder  tube,   '1'  can  be  formulated  by 
substituting  Equation  3.50  into  Equation  3.56.     After  some 
trigonometry , 

=1  -^(|)^/^(l+3/4  |)(3sina/2  +  sin3a/2  +  3cosa/2 
^o 

-cos3a/2)+B  (3.57) 

This  form  of  1'  is  actually  similar  to  that  for  flow  around 
the  edge  of  a  flat  plate  [34,35].     As  is  expected,   for  flow 
close  to  the  rim  of  the  feeder  tube,   the  curvature  of  the  tube 
becomes  less  important.     Fluid  flow  around  the  feeder  tube  thus 
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behaves  much  like  that  around  a  flat  plate.  It  is  also  signi- 
ficant to  note  that  if  Equations  3.50  and  3.56  are  substituted 

_i 

into  Equation  3.55,   lo  can  be  related  to  R  as  w°9R  ^,  a  point 
of  singularity  is  found  for  w  at  R=0 . 


For  fluid  flow  problems,   the  determination  of  the  velocity 
and  pressure  fields    is    important.     The  stream  function  ¥ 
has  been  found  and  expressed  in  terms  of  the  function  g. 
Once  this  stream  function  is  obtained  along  the  stream  line, 
the  velocity  field  can  be  determined  by  using  Equation  3.5. 
However,  the  evaluation  of  the  pressure  field  is  more  compli- 
cated, as  will  be  shown  in  this  section. 

The  pressure  field  can  be  related  to  the  vorticity  by 
using  Equations  3.1  and  3.2  and  the  vector  identity 
V^V=V(V«V)-Vx(VxV) .     Since  u  =  VxV, 

Vp  =  pvVxo)  (3.58) 

Vp  can  be  expressed  in  terras  of  scale  factors  which  are  de- 
fined in  Appendix  A, 


Using  those  relations  between  h,  d  and  y,  and  replacing  ^ 
and  ^2  and  lr,=\i  change  Equation  3.59  to 


Velocity  and  Pressure  Fields 


Vp  = 


(3.59) 


Vp  = 


+  (l-yM*  1^  6y] 


(3.60) 
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On  the  other  hand,  taking  curl  of  w  given  by  Equation  3.53 
yields 


Vxw= 


1  3 


d-y')*  3C 


(c'+y^) 


(3.61) 


The  pressure  gradient  can  then  be  derived  by  substituting 
Equations  3,60  and  3.61  into  Equation  3.58, 


9£  ^  -PV 


9x      (dc^)^(x2+c^')  3y 


C3.62) 


8p  _   pv_ 


9y      (d^^)"(l-y')  9x 


Cx^+yn;') 


(3.63) 


Now,  refer  to  Equation  3.15,  D      can  be  derived  to  be 


=  vdCl-y')     E  G^^_^ix)P^^_^(iM) 


n=l 


(3.64) 


where  ^2n-l^^^  been  given  as  Equation  3.25.  Substituting 

Equation  3.64  into  Equation  3.62  and  Equation  3.63  gives 


9x      (ds^)"^(x'^+c^')     n=l  2n-l 


y  c  '(l-y^)  , 


"(2n-l)  p  (  . 
(x  ) 


(3.65) 
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9P  =     PV  d  ' 
(d;  )  n=l 


2xG2^_^(x) 


(3.66) 


These  equations  can  be  used  together  to  find  p.     In  this 
effort,  say  Equation  3.66  is  used  first  to  integrate  with 
respect  to  y  to  give 


pCx.y)  = 


^     Sn-l(-)^o        o  \     .  dy 


(d^^)  n=l 


(x^+y^C;') 


-  2xG2^-l(-)   ^o     ,  f  '  -2  .  dy  F(x) 

(x  +P  ) 


(3.67) 


This  equation  is,  in  turn,  differentiated  with  respect  to  x 
to  give 


(x  ) 


9x      (dc  )  n=l 


[^2n-l 


(X) 


+2xG 


•  r  ^^  /  ^2n-l^^)  ^ 
2n-l(^)     ^o     ,  2,  .  -S2 

J        (x  +y  ) 


.   8x\^_^Cx)  ,   3  dy 

(x  +V1^Cq  ) 


+  F  (x) 


(3.68) 


Equation  3.68  can  now  be  equated  to  Equation  3.65  to  find 
f' Cx)  at  y=0  to  be 


F  (x)=0 


(3.69) 


which  implies  that 


F( x)=constant 


(3.70) 
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Physically,  this  equation  signifies  that  the  pressure  at  the 
exit  of  the  end  cap  (i.e.,   1J=0)  is  a  constant. 

Equation  3.67  can  be  used  to  find  the  pressure  distri- 
bution.    Note  that  in  this  equation  Gg^  -^(.x)  is  related  to 
g2^_^(x)  by  Equation  3.25,  while  g2n-l^'^^        ^  function  of 
i;^  and  the  Reynolds  number  (through  the  boundary  conditions 
of  g) .     In  view  of  the  definition  of  the  Reynolds  nximber 
(2Vd/v),  the  geometry  of  the  spheroidal  cap  has  thus  been 
totally  accounted  for. 

A  list  of  equations  useful  for  the  integrations  in 
Equation  3.67  is  compiled  in  Appendix  F.  An  integrated 
version  of  this  equation  is  also  given  there. 


CHAPTER  IV 
ANALYSIS  OF  THE  HEAT  TRANSFER  PROBLEM 


Governing  Equation 
Consider  a  two  dimensional  system  in  steady  state.  The 
transport  properties  are  treated  as  constant,   and  the  dissipa- 
tion effect  is  neglected  in  the  analysis.     Under  these  condi- 
tions, the  energy  equation  is  reduced  to 

(V'V)T  =  aV^T  C4.1) 

where  a  =  k/pc  is  the  thermal  diffusivity. 

Making  use  of  the  oblate  spheroidal  coordinates  (x,y)  and 
introducing  the  stream  function  through  the  use  of  velocity 
components  given  in  Equation  3.5  transform  Equation  4 . 1  to 

3T  M  _   9T  91  ^  rj_.     2^  -2     8T-,  ^  A  [(l-uM  ^l"lf4  2) 

3y  3x      8x  9y      ^"^o^ax"-*-^  ^^o       3x 9u  ^  U  y  ; 

which  has  been  derived  in  Appendix  G.     A  change  of      to  4* 
recast  Equation  4 . 2  as 

9T  94'       9T  94-  _  ^o  /  9    .     2,-2     9T-,    .     9  2s  9Tn\ 

-  -^d^  -       \^  )  95Z^      9Tr  ^^^-^  >  d^^j  ('^•2) 

where  Pr=v/a  is  the  Prandtl  number.     Equation  4.3  is  the 
energy  equation  to  be  solved  for  the  fluid  flow  in  the  end  cap. 
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Boundary  Conditions 
For  flows  of  small  Reynolds  number  with  Prandtl  niunbers 
close  to  one,  the  conduction  effect  dominates.     Under  these 
conditions,  the  downstream  condition  usually  has  only  slight 
effect  on  the  temperature  field  in  the  upstream  direction.  A 
set  of  hypothetical  conditions  is  thus  chosen  as  follows 

T(x=0,y)=T^  (4.4) 
T(x=l,vi)=Tg  (4.5) 
T(x,iJ=l)^«>  (4.6) 

Development  of  a  Solution  Technique 
In  an  effort  to  simplify  the  solution  for  the  problem 
given,  a  dimensionless  temperature  is  introduced  as  follows 

T-T 

s 


T  -T 
w  s 


Equation  4.3  is  then  changed  to  be 


y  13^  -       an  =  p7    [-93^  •^(^  -"^0  ^  33^^     9n  '^(^-^  )9]I^)  (^-^^ 


and  for  the  boundary  conditions, 

e(x=0,y)  =  l;   e(x=l,]a)  =  0  (4.8) 

e(x,  =i)^oo;  e(x,p=o)  =  e^(x)  (4.9) 


where    J<y^)  is  the  solution  of  the  one  dimensional  probl 


em , 
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Notice  that  Equation  4.7  is  linear  and  homogeneous  in  6; 
so,  the  sum  of  solutions  to  Equation  4.7  is  also  an  acceptable 
solution.     This  permits  the  use  of  a  perturbation  method  in 
which  the  temperature  6  is  expanded  in  a  parameter  perturba- 
tion as 

00 

e  =  I     (Pr)-^   e  .(x,y)         •  (4. 10) 

j=o  J 

The  unperturbed  part  of  the  solution  is  that  for  j=0,   the  solu- 
tion to  a  pure  conduction  problem.     Equation  4.10  is  a  valid 
solution  to  Equation  4.7  provided  that  the  parameter  Pr  is 
small . 

Substituting  Equation  4.10  into  Equation  4.7  and  equating 
like  powers  of  Pr  results  in  a  differential  equation  that 
relates  subsequent  order  solutions  as 


-r^^^Ji-^-k^'i-i-  J=i.2.3...  (4.11) 

where 


The  boundary  conditions  for  subsequent  9 .   can  be  derived 
in  a  similar  manner  as 


e^(x=0,y)=l,  e^(x=0,y)  =  e2(x=0,y)  =  ...  =0  (4.12) 
e^(x=l,y)=  e-^(x=l,y)  =  ...  =  0  (4.13) 
e^(x,y=l)7^co   ,    e-^(x,y=l)7^  00  .  .  .  (4.14) 
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Zeroth  order  solution 

The  zeroth  order  solution  is  derived  for  a  pure  conduction 
problem  which  satisfies 


8  6  9  6 

and 

86 

6^(x=0,y)  =  l,   e^(x=l,y)=0,    e^(x,y=l)^oo  ,  __^(x,y=0)  =  0 
Solution  to  Equation  4.15  has  been  given  in  [36]  to  be 

oo 

^o  =  n'o  K  Pn(^o-)  ^  \  Pn(^)  (^'l^) 


where        and        are  constants,   and        are  Legendre  polynomials 

p     and  q  are 
^n  ^n 


p^(C^x)  =  i-%(iC^x) 


Equation  4.13  can  be  used  to  establish 


(4.17) 


^n  =  Pn(^o)/^n(^o>  (^-l^) 


It  follows  that 
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The  second  condition  at  x=0  (Equation  4.12)  is  used  to  derive 

A  as 
n 

and        =  A2  =   . . .   =  0  • 
In  this  manner  the  solution  becomes 

(4.21) 


which  can  be  rewritten  as  [36] 


cot~''"(  C  x) 


=  [  1  _^   °     ]/(i  ^^n— )  (4.22) 


cot     C  2cot  C 

o  o 


First  order  solution 

The  first  order  solution  must  satisfy 


with  the  boundary  conditions 


e  (x=0,y)  =  0,   e,  (.x=l,y)=0,   0,  (x,  y=l)^<»  and  e,(x,y=0)  =  0 


44 


Clearly,  the  solution  to  Equation  4.23  now  depends  on  *F 
and  9^,  which  (Equations  3.14  and  4.22)  are  thus  substituted 
into  Equation  4.23  to  give 

'        ^  77~;^r  ;    \   '^'""'^   ,   2,  -S  P2n-l(l^)  ^^'^^^ 

(2cot  n=l  (x  +Cq  ) 

The  left  hand  side  of  Equation  4.24  can  be  expanded  by 
using  the  general  formula  given  in  the  appendix  (Equation  G.3). 
The  resulting  equation  is  multiplied  by  Pj^(y)  a-nd  is  integrated 
from  ]i=0  to  1  to  give 


4?;'  N  ^2n-l^^^  ^ 

  E     n(2n-l)  ,      /    P         (y)P(p)dy  (4.25) 

(2cot  -^C  -tt)  n=l  (x^+C~  ) 


where 


^£  =       eP^(y)dn  •■  (4.26) 


Equation  4.26  is  known  as  the  Legendre  transform  [37].     It  can 
be  noted  that  the  boundary  conditions  chosen  for  d are  all 
homogeneous  and  the  nonhomogenity  for  9  is  all  left  to  the 
9^  term,  see  Table  2.     On  the  other  hand,  the  gradient  of  9^^ 
at  y=0  is  not  known;   yet  this  lack  of  condition  can  be  resolved 
by  changing  the  subscript  I  (for  9  and  P  in  Equation  4.25)  to 
25,-1.     Since  ^25,-1^'^'^  zero,  Equation  4.25  is  reduced  to 


I 
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Table  2.     A  summary  of  boundary  conditions  for   9(x,  y) 


zeroth  order  solution 

first  order  solution 

V=^0     =  0 
o 

0^(x=O,y)=l 

0-,^(x=O,y)=O 

e^(x=i,u)=o 

0^(x=l,y)=O 

96 

0^(x,y=O)=O 

e^(x,y=l)^oo 

0^(x,y=l)?'<» 

total  solution 


e(.x,y)E0^(x,y)  +  0^(x,y) 

satisfies : 

v20=v20  +v20,  =  T-  (  1^      -  1^      )e  (X) 

o         1  ^   9x  9y       3y  8x  '  o^  ' 

EC's: 

e(O,y)=0^(O,y)  +  0^(0, y)  =1+0=1 

0(l.y)  =  0^(l,y)  +  0-1^(1, y)  =0  +  0  =  0 

0(x,O)  =  0^(x,O)  +  0-^(x,O)  =  e^(x)  +  0  =  0^(x) 

0(x,l)  =  0^(7ioo)     +     Q^ifco)     f  CO 
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dG 


where 

(2cot     Cq-^)     n=l  (x  ) 

(4.28) 

Equation  4.28  can  be  further  simplified  by  making  use  of  the 
orthogonality  properties  of  the  Legendre  polynomials;  thus, 

R„     ^  =   '\   !2n^ 

(2cot  -^Cq-ti)     (4n-l)     i^'^+Z^  ) 

Equation  4.27  can  be  solved  by  the  method  of  variation  of 
parameters.     In  this  method,   the  solution  to  Equation  4.27  is 
assumed  to  be 


2n-l  2n-l^    ^  *^2n-l^^o  '         2n-l         ^2n-l     o  ^ 

(4.30) 


where  subscript  £  has  been  changed  to  n.  Equation  4.30  is 
differentiated  to  get 


de2^_,(x) 


Ao„_n(x)  p'      -,(c  X)   +  -,(x)  qA  _-,(C  X) 


dx  "2n-l^^^  ^2n-l^''o"^       "2n-l^^  ^2n-l'''o 

C4.31) 

provided  that 

^2n-l(^)  P2n-l(^o^)       ^2n-l(^>  ^2n-l(^o^)  =  ^  (4. 32) 
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In  the  above  two  equations,  primes  are  used  to  designate 
differentiation  of  the  primed  quantity  with  respect  to  x. 

It  is  possible  to  generate  another  relation  between  A 
and  B  by  substituting  Equation  4,31  into  Equation  4,27,  thus 
giving 


^2n-l(''^ 

P-  x)A'„  -.(x)  +  q'     -.(c  x)B'     -.(x)  - 


2n-l'^f.^  2n-l'^       ^2n-l"  ^o^  2n-l'^  2^~t 
o  X  +c 

o 


(4,33) 


Equations  4,32  and  4.33  can  now  be  solved  as  a  set  of  simul- 
taneous equations  for  A'   and  B'   to  obtain 


A'     ,(x)  =  C^Ro     i(^)qo  x)  (4.34) 

2n-l     ^        o  2n-l^   ^^2n-l  o 

Bl     ,(x)  =  -i;2R„     t(x)p^     ,(5  x)  (4,35) 
2n-l  o  2n-l^  '^2n-l^^o 


These  equations  will  again  be  solved  simultaneously  for  A  and 
B  as  follows. 

Say,  use  is  made  of  Equation  4.35  first,  which  is  inte- 
grated from  x=0  to  x  to  get 


B2n-l('^)  -  B2n-l(0)  =  "^o' (r^2n-l(  ^)P2n-l(^o^)         (4, 36) 


Equation  4.36  can  be  simplified  by  noting  the  fact  that 


e^(x=0,y)=0 


is  the  boundary  condition  for  9^  at  x=0 .  Applying  Legendre 
transform  of  this  condition  provides 
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e2n_i(0)=0  (4.37) 


Since  P2n-1^*^^~'^ '  Equation  4.30  can  be  used  to  get 


B2„_^(0)=0  (4.38) 


Using  this  condition  for  B  reduces  Equation  4.36  to 


^2n-l(^)  =  -^o'^^2n-l(^>P2n-l(^o^)^^  ('^•39) 


In  a  similar  manner  it  can  be  shown  that,   at  x=l, 

0  =  A2n-l(l)P2n-l(^o)  ^  22n-l(  l)^2n-l(  ^o^  ('^•^O) 
Hence,  ^2n-l^^^ 


^2n-l^^o^ 

A2n-l(l)  =  -  pf  '(C  )  ^2n-l(^)  (4.41) 
^2n-l^o 


Use  of  Equation  4.35  is  now  complete;   attention  is  now 
directed  at  Equation  4.34  which  is  integrated  to  be 


A2n-l(^)  -  A2n-l(^)  =  ^o'  ^  R2n-l(^) ^2n-l( ^o^) (4.42) 


The  ^2n-l^^^  °"  left  hand  side  of  this  equation  has  been 

found  earlier  in  Equation  4.41.  Hence,  Equation  4.42  can  be 
rewritten  to  be 


(4 .43) 
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Equations  4.39  and  4.43  are  to  be  used  in  Equation  4.31  to 
f ind  'e2^_^(x)  ,  which  is,   in  turn,  operated  by  an  inverse 
Legendre  transform  to  get 

e  (x,y)  =  Z  (4n-l)[A2^_l(x)P2^_l(C^x)  +  B2^_-L(x)q2^_l(  CqX)]P2^_^(  y) 
n 

(4 .44) 

Equation  4.44  is  the  desired  first  order  solution  to  the 
heat  transfer  problem. 

The  solution  method  presented  above  can  be  repeated  for 
higher  order  solutions.     However,  because  of  the  great  com- 
plexity involved  in  obtaining  the  first  order  solution,  those 
higher  order  solutions  were  not  attempted  in  this  work.  Instead, 
those  constants  A  and  B  in  Equation  4.44  are  re-examined  to 
make  the  first  order  solution  complete. 

Equations  4.39  and  4.43  can  be  expressed  in  terms  of  g 
functions  with  the  substitution  of  Equation  4.29  as  follows 

^2n  l(-)  =   =^1   BC2n^/-%l^  ^Sn-l^^o^^^^ 

2"--^  (2cot-^^-^Tr)   (4n-l)     °  ) 

°  °  (4.45) 


and 

^2n-l^^o^  . 

A2n-l(-)  =  -    ^2n-l(l)  ^   ^1   ""^^^ 

P2n-l(^o)  ^2''°*  ^o-''^  ^^""-^^ 


x  ^2n-l*^''^ 

/        9       2  i(?n^)dx  (4.46) 

1    (x^+c"  )     ^""-^  ° 
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The  appearance  of  B„     ,(1)  on  the  right  hand  side  of  Equation 

Zn- i 

4.46  implies  that  ^2n-l^^^  must  be  found  first  before  ^2n-l^^^ 
can  be  evaluated.     Those  g  functions  in  the  integrands  have  been 
found  earlier  for  the  fluid  flow  problem,  where  as  the  p  and  q 
functions  can  be  expanded  in  series  [36].     The  integrals  in 
Equation  4.45  and  4.46  can  thus  be  carried  out  with  a  numerical 
integration . 


CHAPTER  V 
RESULTS  AND  DISCUSSION 


Data  for  g  functions  were  computed  using  Equation  3.47 
as  shown  in  Figure  7,  all  g  functions  are  of  the  same  order. 
They  range  from  0.5  to  2.5.     Recall  that  the  magnitudes  and 
slope  of  the  g  functions  at  x=0  are  governed  by  the  boundary 
conditions  listed  in  Table  1,     Since  these  are  curves  for  the 
same  Re  number,  the  similarity  of  the  g  curves  of  Figure  7 
is  not  unexpected. 

Recall  that  the  conditions  for  g^  and  g^^  at  x=0  have 
been  established  in  Chapter  III,  as  g^  -j^=Cl/9  )g-|^  -j^ ,  and 
inferrence  was  drawn  at  that  time  that  the  g  functions  should 
converge  rapidly.     This  point  is  now  substantiated  by  the 

curves  shown  in  Figure  7.     In  fact,   gcCx)  for  i;  =2.5  is  found 

o  o 

to  be  everywhere  zero.     Finally,  there  appears  to  be  little 
change  in  the  trends  evidenced  by  the  g  functions  when  the 
Re  number  is  changed  from  0.1  to  0.3,  see  Figure  8. 

Stream  Functions 
The  g  functions  permit  the  evaluation  of  the  stream 
function.   Equation  3,15.     A  plot  of       for  Re= .  1  and  i;^=0.5,  2, 
and  3  is  presented  in  Figure  9.     Notice  that  Re  has  been  de- 
fined in  this  work  as  2Vd/v,   and  5^=c/d.     d  is  the  radius  of 
the  feeder  tube,  whereas  c  is  the  depth  of  the  cap.     A  large 
thus  refers  to  an  end  cap  having  a  large  annular  passage 
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g(x) 


Figure  7.     S2n-1  ^^^"^^^  Re=0 . 1 
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relative  to  the  feeder  tube.     Such  a  condition  is  conducive 
to  a  circulation  flow  near  the  rim  of  the  tube,   and  the  stream 
lines  in  Figure  9Cc)  and  10  substantiate  this  point.     A  magni- 
fication of  the  flow  pattern  near  the  rim  of  the  feeder  tube 
shown  in  Figure  lOCb)  is  provided  in  Figure  11. 

Velocity  Fields 

The  velocity  field  was  determined  using  Equation  3.5. 

The  velocity  distribution  was  plotted  in  three  dimensions 

and  shown  in  two  views  for  c  =2  and  Re=0.1  in  Figure  12.  At 

o  " 

first  sight,  the  entrance  velocity  in  the  figure  may  appear 
to  be  high  compared  with  that  at  the  exit .     This  turns  out 
to  be  a  matter  of  geometry.     For  the  end  cap  under  investiga- 
tion in  this  plot,   C^=2=Cc/d)  is  equivalent  to  an  (a/d)  ratio 
of  2.236;  the  exit  velocity  from  the  end  cap  is  thus  very  small 
in  order  to  satisfy  continuity  of  flow.     Also  notice  that  this 
Cq=2  is  close  to  the  limit  beyond  which  a  circulation  will 
appear  near  the  rim  of  the  tube. 

Pressure  Fields 
The  pressure  drop  in  the  end  cap  is  expected  to  be  large 
because  the  flow  is  deflected  180  degrees  inside  the  cap.  It 
is  informative  to  examine  how  the  pressure  distributes  in  the 
flow  region. 

The  pressure  field  was  determined  using  Equation  3.67  and 
plotted  in  three  dimensions  for  5^=2  and  Re=0.1,   see  Figure  13. 
The  pressure  drop  in  found  to  be  quite  marked.     There  is  a  large 
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Figure  11.     A  Magnification  of  the  Streamlines  Near  the 
Edge  of  the  Feeder  Tube  for  r,^=3  and  Re=0 . 3 
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x=l  pc^/pv^d=1.53 


Figure  13.     Two  Views  of  the  Pressure  Distribution 
in  the  End  Cap  for  C  =2  and  Re=0 . 1 
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pressure  exerted  on  the  surface  of  the  end  cap  in  the  axial 

direction.     This  is  certainly  a  result  of  the  conversion  of 

the  velocity  head  to  pressure  head  as  the  axial  velocity 

decelerates  to  zero  when  the  flow  impinges  on  the  cap  surface. 

The  velocity  drop  in  the  axial  direction  has  also  been  computed 

2 

and  tabulated  as  uc  /vd  in  Table  3. 

Temperature  Fields 
Equation  4.10  was  used  to  plot  the  temperature  field  in 
Figure  14.     Conditions  chosen  for  this  plot  are  ?q=2  and  Re=0.1, 
The  solid  line  in  the  figure  represents  the  pure  conduction 
case,  that  is,  the  zeroth  order  solution  of  the  temperature 
distribution  in  the  perturbation  method.     An  addition  of  con- 
vection effect  displaces  the  isotherms  to  the  right,  which  is 
not  unexpected.     Notice  that  for  the  problem  solved,  the 
temperature  is  truncated  to  the  first  order;  yet  the  conduc- 
tion effect  is  dominant  as  is  clearly  seen  in  the  plot. 

Local  Convective  Coefficient 


At  the  surface  of  the  cap,   the  heat  transfer  coefficient 

=0.     In  terms  of  the  Nusselt 


is  defined  as  h(T  -T  )+k 

w     s  9n 


s 

number  Nu  and  the  dimensionless  temperature  G,  the  above  defi- 

9  6 

nition  becomes,  Nu  = 

x=l 

Equations  4.22  and  4.44  were  used  to  determine  the  Nu(y) 
plotted  in  Figure  15  for  Pr=0.5,   0.6,   0.8,   0.9.     The  conditions 
are  given  in  the  legend.     Notice  that  the  Nusselt  number  de- 
creases with      and  approaches  the  limiting  case  of  a  pure 
conductive  situation  of  Nu=0.362.     The  Nusselt  number  is  high 
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Figure  15.     Nu(y)  Curves  for  Cq=2  and  Re=0,l 
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in  the  region  where  the  wall  pressure  is  high.  This  effect 
becomes  less  prominent  if  the  Prandtl  number  is  low. 


CHAPTER  VI 
CONCLUSIONS  AND  RECOMMENDATIONS 


Conclusions 

Based  on  the  work  done,  one  can  draw  the  following 
conclusions:  . 

1.  An  exact  solution  has  been  developed  for  fluid  flow 
and  heat  transfer  in  the  end  cap  of  a  single-pass,  return  flow 
heat  exchanger  for  the  special  case  of  steady,   low  speed  flow 
condition.     In  this  case  the  product  of  velocity  components 
can  be  neglected  in  the  analysis.     In  the  solution,  the 

stream  function  is  expressed  in  a  series  of  Legendre  polynomials 
in  which  only  three  terms  are  found  to  be  significant .  The 
series  converges  rapidly  for  the  geometry  and  the  Re  number 
range  considered  in  this  study. 

2.  A  parameter  perturbation  method  is  used  to  solve 
for  the  temperature  filed  in  the  end  cap.     The  zeroth  order 
perturbation  solution  corresponds  to  a  pure  conduction  problem. 
Because  of  the  great  complexity  involved  in  the  solution,  the 
perturbation  is  developed  only  up  to  the  first  order  with  a 
hypothetical  case  of  uniform  inlet  and  wall  temperature  condi- 
tions as  boundary  conditions. 

3.  In  the  fluid  flow  study,  the  obtained  stream  functions 
are  used  to  determine  the  velocity  and  pressure  fields.  A 
recirculation  flow  is  found  near  the  rim  of  the  feeder  tube 
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when  the  tube  annulus  region  has  a  cross  section  that  is 
large  compared  to  that  of  the  tube. 

4.     For  the  heat  transfer  study,  the  conduction  effect 
is  found  to  be  dominant.     Larger  convective  coefficients  are 
found  in  the  wall  region  where  the  pressure  is  high. 

Recommendations 
Some  recommendations  are  made  as  follows: 

1.  It  is  clear  that,  because  of  the  great  complexity 
involved  in  the  exact  solution,   only  some  simplified  cases 
can  be  studied.     For  those  cases  whose  solutions  are  given 
in  this  work,   further  numerical  solutions  or  experiments  are 
needed  for  verification. 

2.  Whether  a  higher  order  perturbation  is  needed  in  the 
solution  of  the  heat  transfer  problem  is  also  subjected  to  a 
rigorous  verification. 

3.  The  range  of  the  parameters  (.Re  and  Pr)  studied  in 
this  work  is  quite  limited.     It  is  necessary  to  study  their 
limits  under  which  this  analytical  solution  is  useful.  This 
can  only  be  made  by  comparing  the  results  obtained  by  numerical 
solutions  of  the  problem. 


APPENDICES 


APPENDIX  A 
DERIVATION  OF  EQUATION  C3.8) 

It  is  convenient  to  establish  relations  between  the 
cylindrical  coordinates  (r,z,({))  and  the  oblate  spheroidal 
coordinates  (n,6,(j))  with  the  help  of  the  cartesian  coordi- 
nates (x,y,z).     In  terms  of  the  cylindrical  coordinates, 
the  cartesian  coordinates  can  be  expressed  as  [38] 

where,   ^i"^^'   C2=cos(})  and  ^3=2- 

In  terras  of  the  oblate  spheroidal  coordinates  (r\,Q,(^),  the 
cartesian  coordinates  become 

X  =  ^3\/(^i+d^)(l-^2) 

y  =\/(q+d2)(l-^Kl-^) 

where,   5,=dsinhri,    ^,,=0036  and  ^o=cos4). 
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The  relations  between  the  cylindrical  and  the  spheroidal 
coordinates  can  then  be  established  by  relating  these  two 
sets  of  x,y  and  z  equations.     It  follows  that 

r=d  coshri  sinG 
z=d  sinhri  cosO 

By  defining,   cos9=y  and  sinhn=  C ,   it  is  possible  to  rewrite 
r  and  z  as 


which  relate  (r,z)  to 

The  introduction  of  the  curvilinear  coordinates  developed 
above  enables  a  speedy  solution  to  the  fluid  flow  problem. 
Such  a  transformation  is  not  uncommon  in  the  field  of  mathe- 
matical physics  in  the  solution  of  punch,   crack  and  torsion 
problems,  among  others  (e.g.,  [39]). 

It  is  easy  to  show  that  the  surfaces  y^constant  and 
5=constant  are,   respectively,   hyperboloids  of  revolution  and 
oblate  spheroids.     The  spheroid  r;=0  is  the  circular  disk  r<d 
and  z=0,  while  y=0  gives  the  circular  ring  surrounding  the 
disk  (Figure  A. 1)  . 

In  fluid  mechanics  the  stream  function      satisfies  the 
continuity  equation  which  can  be  expressed  in  vector  notation 
as  follows 


r=d(l+C^)^(l-y^)^ 


(A.l) 


z=d;y 


(  A.  2) 


V  •  v=o 


(A. 3) 
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Figure  A.l.     Oblate  Spheroidal  Coordinate  System 
  :  Fluid  Flow  Region 
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Use  of  axisymmetry ,  the  divergence  of  V  may  be  written  as 

'  •  ^  =  h^l^  ^(hghgV,)  .  3f^(h3h^V,)]  (A.4) 


where  h    are  scale  factors,  defined  as 
a 


9  ^^i  2 

1  a 


where  x.   are  the  cartesian  coordinates. 
1 

These  factors  take  the  following  forms  in  the  oblate  spheroidal 
coordinates : 


h 


1   \i  q+d^ 


^2 


^3 


'(C^d2)(l-5') 


Substituting  h^  into  Equation  A.4  and  making  use  of  Equation 
A. 3  gives 


^  lai+d^q)^a{+d^)h^}  +  ^  [(q+d2^)*(i-cp*v2]  =0 

Since  i-^=dc,  and  ?2~^'  above  equation  can  be  recast  in 

terms  of  (C,y)  coordinates  as 
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It  is  now  possible  to  define  the  stream  function  i) , 
which  is  related  to  the  velocity  components  u  and  v  as 
follows 


u  =   ^1  J.  ^  (A. 5) 


v  =   I  (A. 6) 

d^;^+U^)^(l-y^)^ 


Clearly,  the  continuity  equation  is  satisfied  by  ^  identically. 
These  velocity  components  have  been  used  to  write  the  velocity 
V  in  Equation  3.5. 

The  derivation  given  in  Chapter  III  for  the  momentum 
equation  involves  use  of  the  definition  of  vorticity  given  as 

0)  =  VxV 


In  terms  of  the  curvilinear  coordinates  ^3^ '  "^^^^ 

becomes 


In  the  (?,y)  two  dimensional  coordinate  system  that  is  of 
interest  in  this  work, 


^  =  (i-y^)^(i+c^)^  ;  J_ 


(i-y^) 


V 


_8_ 

8y 


u 


2  \  2 


(1+r ) 
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This  w  equation  can  now  be  expressed  in  terms  of  the  stream 
function  as 


0) 


1  ^     1  all 


(A. 7) 


Finally,   in  order  to  obtain  an  expression  for  the  govern- 
ing equation  in  (C,y)  coordinates,  the  curl  of  Equation  A. 7 
is  repeated  twice,  which  gives 


d^  Uc'+y') 


(l,,a)  all] 


(1-y^)  9' 


i2       9y2  + 


9C' 


9U' 


(A. 8) 


A  close  examination  of  this  equation  suggests  the  defini- 
tion of  an  operator 


9y' 


which  permits  a  rewriting  of  Equation  A. 8  as 


[   D^^]  =0 


(A. 9) 


An  alternative  expression  to  Equation  A. 9  is  to  consider 
use  of  dimensionless  quantities  'I'  and  x. 


'      vd  '  - 
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Equation  A. 9  is  thus  changed  to 


T)^[—-^          D^^]=0  (A. 10) 

X  +y 


where , 


APPENDIX  B 

DERIVATION  OF  EQUATIONS  3.21  THROUGH  3.23 

Let  f(y)  in  Equation  3.13  be  a  function  that  satisfies 
Dirichlet  conditions  in  the  interval  (-1,1).     Then,   f(y)  can 
be  expanded  as 


f(y)  =  (1-y^)     E     B2^_,  Psn-i^^^  (^-^^ 
n=l 


Invoking  use  of  the  orthogonality  characteristics  of 
(1-y^  )P2j^_-,^(y  )  over  C-1,1),   it  can  be  shown  that  [33] 


In  the  special  case  of  a  fully  developed  flow  at  the 
exit  of  the  feed  pipe,  the  velocity  is  related  to  the  radius 
according  to  [3] 


=  u=  2V[1  -  C^)^]  (B.3) 


Since  at  the  exit  plane     5=0,  Equation  3.3  can  be  used  to 
rewrite  the  velocity  u  as  a  function  of  y,  that  is, 

u  =  2V[l-(l-y^)]=2Vy2  (B.4) 
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In  addition,  u  is  related  to  i|j  according  to  Equation  A. 5. 
At  the  location  C=0, 


u=-^^  (B.5) 
d  y  8y 


Equations  B.4  and  B.5  can  now  be  combined  to  give 


•4"  =  2Vy2  (B.6) 
d  y  9y 


and  ijj(0,y)  derived  as 


iP(0,y)  =  -V  d^  ^  +  c  (B.7) 


The  constant  c  of  integration  can  be  determined  from  the 
fact  that,  along  the  axis  of  symmetry  (i.e.,  y=l),  the  stream 
function  is  zero.  Hence, 


_  2 

iKO.y)  =  ^  U-y"*)  (B.8) 


In  terms  of  Reynolds  number  Re=2Vd/v,  this  ij;  can  be  changed 
to  as 

H'(0,y)=Re(l-y'*)/4  (B.9) 

At  this  stage  it  is  possible  to  use  the  above  equation 
in  Equation  B.l  to  determine  62n-l'     ^^^^  '^^^  help  of  Equation 
B.2,  it  can  be  deduced  that 


77 


3^  =  3/10  Re 
63  =  1/30  Re 

and 

3^      =0  for  n  >  3  (B.IO) 

2n-l  — 


Substituting  these  in  Equation  B.l  gives 

'F(0,y)  =  fCy)=4Re(l-y2)[3/5  p|(y)  +  l/15  PgCy)]  (B.ll) 

Comparing  Equation  B.ll  with  Equation  3.15  enables  the  g 
function  to  be  found  immediately  as 


g^(0)  =  3/10  Re 
ggCO)  =  1/30  Re 

and 

^2n-l(0)=0  forn>3  (B.12) 


APPENDIX  C 
DERIVATION  OF  EQUATION  3.36 

The  governing  equation  to  be  solved  for  the  fluid  flow 
problem  is 


^     ,  D^'i'  =  0  (C.l) 


where , 


A  solution  ^  has  been  suggested  in  the  form 
n=l 

By  using  the  identity 

[(l-y2)p^^_^(3j)j"  +  2nC2n-l)P2n_iCp)  =  0  (C.3) 
and  defining 

G2n-l(^)  =  (^^^^o'^S2n-l(^^  "  2n( 2n-l )g2^_^( x)  (C.4) 

it  can  be  derived  that 

D^v}/  =   (l_y2)      ^     G2^_;^(x)P2^_^(y)  (C.5) 
n=l 
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Equation  3.24  can  be  derived  by  subsituting  Equation  C.5 
into  Equation  C.l  and,   after  some  manipulation,  this  gives 


n=l 

-  [4^;%[(l-yMP2^_^(y)]'  +  2[(n+l)(2n-3)(x2+yH^V6yH^'-2c;'] 


(X)  =  0 


(C.6) 


In  order  to  sort  out  G  from  this  equation,  the  identities 

3.26  and  3.27  are  substituted  for  those  underlined  terms  in 

Equation  C.6,  and  after  some  algebra,  Equation  C.6  is  trans- 
formed to 


n=l 


2     -2  8n'^-4n-3 
^  ^o  (4n+l)(4n-3) 


^2n-l  -4^(^^-^^o')4n-l 


.  /  8n(2n-l)-2(2n^-n+3)(8n^-4n-3)  _ 
"^V  (4n+l)(4n-3)  ^o 


(n+l)(2n-3)x^-2^  ^ 


V  r.        /   N  ^        2^^-2.       2n(2n-l)        -"^  p„ 
^2n-l?  P2n-l(^)  ""^o  >  (4n-l)(4n-3)  ^o  ^2n-l 


4n(2n-l)(2n+l(n+l)  - 
(4n-l)(4n-3) 


o 


Sn-1?  ^2n-3(y) 


2^^-2.  2n(2n-l) 


_2 

G, 


o  '   (4n-l)(4n+l)  ""o  2n-l 


4nC2n-l(n-l)(2n-3)  -2 
(4n-l)(4n+l)  ^o 


Sn-l?P2n+l(^)  =  0 


(C.7) 


A  general  equation  for  G  can  then  be  derived  from 
Equation  C.7  by  (i)  expanding  the  series  out,   (ii)  grouping 
terms  of  p|(y),  PgCy ) ,    ...    ,   and  finally  (iii)   equating  each 
coefficient  to  zero.     The  differential  equation  for  G  is  given 
as  Equation  3.28. 
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Following  the  procedure  described  in  the  text  to  derive 
Equation  3.36,   it  can  be  shown  that  the  coefficients  A  to 
Q  in  this  equation  can  be  related  to  a,  b,   and  c  (Equations 
3.29  through  3.35)  by 


(C.IO) 


D,„_,=-2„(2n-l)b(l)^  (C.ll) 

E2„.,=-2„(2„-l)b(°\  (C.12) 

H2n-r<'''*^o'>4n^l 

■'2n-r-2(n-2)(2n-l)a<2)3  .  (x^-H?;^  )a(0>,  (C.15) 
«2„-r''-2n'l 

Q2„.i=-2(n+l)(,2n+l)4°|j  (C.20) 


APPENDIX  D 

SOLUTION  OF  EQUATION  3.36  BY  A  FINITE  DIFFERENCE  APPROACH 


The  set  of  equations  to  be  solved  for  g  is 


*2n-l4n-l  *  ■^2n-lS2A-l  *  ""-^n-A-l  *  °2n-lS2n-l  *  '^2„-l'^2n-l 


^2n-l""2n-l4n-l       ^2n-1^2A-3  ^  ^^211-1^211-3  ^2n-1^2n-3 

for  n=2,3,    ...    ,  N  (D.2) 

S2n-rL2n-l4n  +  l  ^  ^^2n-lS2A+l      ^s^.^g^^^^  +  ^211-1^211  +  1 

for  n=l,2,    ...    ,  N-1  (D.3) 


2n-l 


Coefficients  A  through  Q  in  these  equations  have  been  given 
in  Appendix  C. 

Using  the  central  difference  scheme  given  in  [40]  permits 
writing  the  derivatives  of  g  as 


g'(xj)=  ^(-g(x._^)+g(x.^^)) 
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g"(Xj)=  ;^(g(x^_^)-2g(x^)+g(Xj^-^)) 

g"'(xj)  =  -^(-g(x^_2)  +  2g(Xj_^)-2g(Xj^^)+g(x.^2)) 
g^^(Xj)  =  ;^(g(x._2)-4g(Xj_^)+6g(Xj)-4g(Xj^-^)+g(x_.^2>) 

where  x^  denotes  position  at  index  j,   and  h  represents  the 
step  size  in  x-direction.  Thus, 

Xj=(j-l)h 

2 

The  derivatives  of  g  given  above  are  accurate  to  h  according 
to  Taylor's  expansion.  The  relation  between  x  and  j  is  shown 
in  Figure  D . 1 . 

Equations  D.l  through  D.3  will  be  used  to  write  g  rela- 
tions at  point  X..     It  is  thus  convenient  to  use  double  sub- 

3 

scripts,  the  first  one  identifying  the  term  in  the  series 
expansion  in  Equation  3.15,  the  second  one,  the  position  of 
the  point,  for  example, 

^2n-l^^j''"^2n-l,  j 

Substituting  Equation  D.4  into  Equation  D.l  gives 

(A2n-l,J-^2n-l,j  |)S2n-l .  j-2-^( -^A2n-1 ,  j-^^2n-l ,  j  ^^Sn-l.j 

3 

-°2„-l,J  V>82n-l,J-l*<6A2„.^_  .-2C2„.i  J  h=^+E2n-l,J  """^Hn-Li 
*(-«2n-l,r^2„-l,J  h-Sn-l.j  h'+Dgn-l.j  V^2„-l,j.l 
*<*2n-l,J™2n-l,J  2)'^2„-l,J+2  =  -«2n-l.j  ''""^an-LJ  (D.5) 


Figure  D.l.     g(x)  in  the  Neighborhood  of  a  Point 
and  the  Discretization  Scheme 
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Likewise,  Equation  D.2  is  changed  to 

(H2„-l,j-l2n-l,J  l>S2n-3,j-2  ■^<-«2n-l , /^2„-l ,  j  ''*'2n-l  .j"^^ 
82„-3,J-l  +  («2n-l,j-2J2n-l,j  "Xn-lJ  ""jgsn-S.j 
*(-''«2n-l,J-l2„-l,j  ''*''2n-l,j  '''>e2n-3  .        *<"2n-l , /^2n-l ,  j  I' 
S2n-3,J.2  'Hn-l.i 

and  for  Equation  D.3, 

(L2n-lJ-^^2n-l,j  |)^2n-.l,j-2  ^(-4L2n-l ,  j^^Wl ,  J  ^^N2n-l,J 
^2n+l,J-l  -^(^WLJ-^Nsn-lJ  ^Hn-l,J  h^)S2n+l.j 

^(-^L2n-lJ-^2n-l,J  ^  ^  N2n-l,j  ^'^Ssn+l.j  +  l  •^(L2n-1 ,  j^^2n-l .  j  l> 

s  =S  h'^  (D.7) 

^2n+l,j+2  ''2n-l,j 

The  boundary  conditions  g'(0)=0  and  g'(l)=0  can  be 
written  as 

g2n-l,0=S2n-l,2  g2n-l , M-l=S2n-l , M+l  ^^'^^ 

respectively,  and  n=l,2,3,    ...    ,  N. 

By  taking  j  to  be  2  through  M-1  and  making  use  of 
Equation  D.8  it  is  possible  to  write  Equation  D.5  in  a  matrix 
for  the  nth  term  in  the  series  expansion  3.15.     That  is. 
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0  0... 

^21  ^22     ^23  ""^24 

"•^Sl  ^32     ^33  ^34 

°  ''"m-2  M-4'  •  •  •  •  • 

"""m-I  M-3-  •  • ' 
0       0...  0 


35 


0 
0 

0..  0 

"'^M-l  M 
0  0 


^2n-l,l 
^2n-l,2 
^2n-l,3 


0 

4  4 
■^n-1,2^  ~^2n-l,2^ 


g. 


'2n-l,M-l 


l^2n-l,M 


■^n-l  ,M-l^  "^2n-l,M-r 
0 

-  (D.9) 


which  is  valid  for  n=l,2,  ...  ,  N.  The  elements  in  the  coef- 
ficient matrix  are  given  as  follows 


T22=^^2n-l,2-^2n-l,2  I  -2Sn-l,2  ^'•^E2n-1,2 

h  2  ^ 

"^M-l  M-r'^^2n-l,M-l''^2n-l,M-l  2^Sn-l,M-l  ^  ■'^2n-l,M-l  ^ 


^33 


21 


='m-2  M-2=6A2n-l,J-2C2n-l,j^'^^2n-lJ^'   '   ^  =  ^ '  ' 


=  Tm-1  M-2  =  -^V-l,j^^2n-l,/-'Sn-l,j^^^°2n-lJ 


M-2 
3 


j=2,         ,  M-1 


23 


••  =  Vl        -^^2n-l,j-^2n-l.j^^Sn-l.J^^^°2n-l,j  V' 


3=2,    . .    ,  M-1 


31 


24 


•••   =  Vl  M-3=A2n-l.J-^2n-l,j   2'     3  =  3,    ....    ,  M-1 


"  "^1^1-2  M  ^2n-l,  j'^^2n-l.  j  2'      ^  ^' 


,  M-2  (D.IO) 
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In  a  like  manner,  D.6  can  be  expressed  as 


0. 


^23       '^24  ° 


'^31     '^32       '^33  ^^34 


'35 


M-2  M-4 
R 


R 


.  .  0 

0 

0.  .  0 

M-2  M 


M-1  M-3     ■  ■    '^M-l  M 


^2n-3,l 


0  ^ 

4 


^2n-l,2^ 


I ^2n-3,M 


^2n-l,M-l^ 


(D.ll) 


for  n=2,3,    ...    ,  N.     The  elements  in  the  coefficient  matrix  are 


f^22=7«2n-l ,  2-l2n-l .  2'^-^'2n-l ,  2h'^^2n-l .  2^^ 

h  2  ^ 

'^M-l,M-r'^^2n-l,M-l"^^2n-l,M-l  2"^'^2n-l,M-l  ^  "^Sn-l.M-l^ 


"^33" 


Rm-2  M-2=6«2n-l,J-2J2n-l,/-^Sn-l./'   '   J  =  3,...M-2 


^^21" 


^M-l  M-2=  -^«2n-l,3"^2n-l,jh^J2n-l,J^     '  J=2,..,M-1 


■^23" 


^M-l  M=  -^«2n-l,j-^2n-l,J^^'^2n-l,j^^ 


3=2,.. ,M-1 


^^31" 


^M-l  M-3=«2n-l,j-^2n-l,J  I  '     ^  =  ^ '    ■    ■  '^^"^ 


^^24" 


'^M-2  M"^2n-1,  j'^^2n-l,  j  2   '   ^  ^' 


,M-2  (D.12) 
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Finally,  Equation  D.7  is  expressed  as 


0   .     0.  . 

^21     ^22     ^23     ^24  ^" 
^31     ^32     ^33     ^34     ^35  ° 


0 
0 
0 


0, 
0, 


^M-2  M-4    ^M-2  M 

^M-1  M-3" '^M-l  M 


0. 


^2n+l , 1 


S2n-1.2^ 


0  ^ 

4 


= 


^2n+l,M 


^2n-l,M-l'' 
0 

(D.13) 


for  n=l,2,    ...    ,N-1.     Elements  S  are 


S22=^Sn-l ,  2-^^2n-l ,  2  h^'^2n-l ,  2^^^Sn-l .  2^^ 

h  2  < 

^M-l,M-r'^^2n-l,M-l'^^2n-l,M-l  2~^^2n-l ,  M-1^  "^Sn-1,M-1^ 


S  = 
^33 


M-2=-«^2„-l,J-2N2„-l,3h'-^«2n-l,j'^^  •   3=3,.., M-2 


21 


=  Vl  M-2=-4L2„-l,/^2n-l,3h^N2„-l.j^     '   J  =  2,...,M-1 


S 

^23 


=  S,,  ,     =-4L„^  ,    .-M„^  ,    -ih+N  h     ;  j=2,...,M-l 

M-1  m        2n-l,j     2n-l,j  2n-i,j 


^31- 


M-3=''2n-l,j-^2n-l.j  I  •  J^a.--.*"-! 


^24= 


■Sm-2  M=4n-l,/"2n-l,J  2   '  J==  ""^ 


(D.14) 
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Those  three  matrix  equations  given  above  can  now  be  combined 
to  give 


0.  . 

0 

V 

-1 

ti]n=2 

0 

1 

V 

-3 

V  • 

-2N-3 

0 

0 

=N^--'n=N 

V 

-2n-l 

=  0 


(D.15) 


where , 


V  = 
-1 


^1,1 
■1,2 


g 


g 


3,1 
3,2 


g 


l.M  J 


g 


3,M 


(D.17) 


and  so  on  .    .  . 

Components  g  given  in  _^  vectors  can  be  further  simplified 
by  using  those  boundary  conditions  which  have  not  been  used  so 
far.  Thus, 


g^(0)=3/10  Re         is  changed  to 


g  (0)=l/30  Re  becomes 


g^  -,^=3/10  Re 
gg  -^=1/30  Re 


CD. 18) 
CD. 19) 


g2^_-|^(0)=0  for  n  ^  3  becomes 

S2n-1^^)^°  becomes  ^1 ,M=^3 ,M= • • =°  ^^'^^^ 

The  computer  program  is  compiled  in  Appendix  H. 


APPENDIX  E 
DERIVATION  OF  EQUATIONS  3.50  AND  3.51 


Equations  A.l  and  A. 2  in  Appendix  A  can  be  used  to  derive 
as  functions  of  r  and  z  as  follows 


,  =     1       Ur2dWz2-d^+  wW+dV:dV^d6z2}  "  (E.l) 
I  =  2*dz  {(r^d^d^z^-d^  -.sJ^^^TTTJ^TT^^^}  * 


2^d 

y  = 


(E.2) 


Consider  cylindrical  coordinates  (r,z)  for  a  moment. 
These  coordinates  can  be  related  to  R  and  a  based  on  the 
geometry  given  in  Figure  6.  Hence, 

r=d+Rsina  (E. 3) 

z=R  cosa  CE.4) 

Substituting  Equations  E.3  and  E.4  into  Equations  E.l  and  E.2 
gives 

1  1 
C  =(f)^  ^  [(R+2d  sina)+(R^+4d^+4d  Rsina)*]  (E.5) 

y  =(2R)*  cosa[(R+2d  sina)+( R^+4d^+4d  Rsina)*]  (E.6) 

The  expressions  inside  the  braces  in  these  two  equations 

are  now  designated  as  A^(R>  (for  C)  and  A^{R)   (for  y),  and 

Taylor's  expansions  are  used  to  write 
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p  4  1  a'(0)        a''(0)  2 

C  =  (f)     ^  [AiCO)  +  R  +  -il—        +  (E.7) 

4  a'(0)        a"(0)  2 

y  =  (2R)     cosa[A2(0)+  -fi        R  +  -|l        R    +   ..J  (E.8) 


For  A    and  A„  to  be  expanded  about  R=0.       These  A  and 
their  derivatives  can  then  be  found  and  the  C  and  y  equations 
rewritten  as 


,  =  (f)*  \  ((2d)*  .  R  .  ^^-^'p^  ^  ^  ..  )(l+sina)* 

^       ^    I  2(2d)^  4(2d)^'*/'^    2  J 


(E.9) 


rnr.^i                                   1           „  ^  (l+2sina)         ,  "1 
y  =  (2R)^  cosa  {(2d)  ^  oj^  R  +  ^  eto^  —       •  •  > 

I  2(2d)'^/^  4(2d)^/^      2  J 


(1+sina)  *  (E.IO) 


With  the  trigonometric  identities 


and 


(1+sina)*  =  sina/2  +  cosa/2 


cosa(l+sina)  *  =  cosa/2  -  sina/2 


Equations  E,9  and  E.IO  can  be  recast  as 

1   /•  2  \ 

C  =  (|)^  (l-2sina)  +  ...  \  Ccosa/2  +  sina/2) 


(E.ll) 


rR^*  ^^_i/-R^4-A.('R^^(l■^2sina)  +  ...  \(  cosa/2  -  sina/2 


(E.12) 


These  are  Equations  3.50  and  3.51  given  in  the  text. 


APPENDIX  F 
DERIVATION  OF  PRESSURE  FIELD  EQUATION 

The  integrals  in  Equation  3.67  are  carried  out  in 
this  appendix.     Equation  3.67  can  be  written  for  F(x)=0  as 

p(x,y)(dc^)-^/(pv^d)=     E     G        (X)  /   -^f-^^-^  dy-2xG  (X) 
o  n=l    ^"  (X  ) 

°  (x^+y^C;')^ 

where  the  series  terminate  at  N=3.     ^2n-l  ^2n-l 
tions  are  given  as  follows: 


^2n-l(^)=(^^^^o')  g2n-l(^>-2"(2n-l)g2^_^(x) 

G2n-l(^)  =  (^^-*-^o')^2A-l(^)^2xg^n-l(^)-2"(2n-l)g^^_^(x) 

Pi(y)=i 

P3(u)=3/2(5y2-l) 

P^(y)  =  15/8(21y'*-14y2  +  l) 


Integrals  in  Equations  F.l  can  be  carried  out  with  the 
help  of  the  formulae  given  in  [41] .     There  is  derived  that 
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p(x,y)(dC  )^/(pv2d)=[-  -|  +  (-^  -  5c^)3/2G3+(63x2c^) 

X  X 

+14C^  -  l/x^)15/8G5  +  G^/x  -         +  5xC^)3/2  G3  +  (21x^C^ 
+14cJx+l/x)15/8  Gg]CQtan~^(y/xC^)+[-G-^/x+(i  +5x^^)3/2  G3 

+  (2lxyn'Q-14C^x-l/x)15/8  G^  ]  y/Cx^+y^c^') 
+{G3+[21(^  x2cJ)-14]G5/4  -  63/4  xC^g^}  15/2y^J  (F. 


APPENDIX  G 
DERIVATIO^J  OF  EQUATION  4.2 


Both  sides  of  the  energy  equation 

(V  •  V)T  =aV^T  (G. 1) 

can  be  expressed  in  terms  of  the  general  curvilinear  coordi- 
nates (K-^,   ?2'  ^3-^ 

V  V 

(V  .  v^T  =  —  +  —  - —  (G.2) 

and 

,           .  h^h„                 r,        h.ho  T 

=  E^i^         C  ^  ^)  .  4(  ^  31^)] 

where        and        are  velocity  components  of  V,   and  h^,  and 
hg  are  those  scale  factors  defined  in  Appendix  A. 

Upon  substituting  of  Equations  G.2  and  G.3  and  making 
use  of  those  h  equations  mentioned  above  changes  Equation  G.l 
to 


(^?+d2)*v,       +  (i-a)*v2  If-  =  ^  .{  ^  la\+d^)  — ] 

^  1  9?_^  2       2  (Cj+d^l)^     95^         ^  9^;^ 

+  [(1-^)  — ]   )  (G.4) 


3^2 
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Notice  that  the  curvilinear  coordinates  (^-j^, 
related  to  (x,y)  by 


and  for  the  velocity  components, 


V-|^=u    and  V2=v 


Equation  G.4  can  then  be  expressed  in  terms  of  these  new 
variables  as 


(X  +C^')^u—  +  (l-y^)^v  —  =  ^        o        -2  1  ^aSE  '-(^  -"^o  )  93EJ 

Upon  using  the  relationships  between  the  velocity  com- 
ponents and  the  stream  function  (Equations  A. 5  and  A.6), 
Equation  G.5  can  be  changed  to 


f  S  -      f  =  «^^„-<^  i^^'<^  [a-.^)  (0.6, 


which  is  Equation  4.2, 


APPENDIX  H 


COMPUTER  PROGRAM  AND  SAMPLE  DATA 
FOR  STREAM  FUNCTIONS, 
VELOCITY  AND  PRESSURE  FIELDS 


97 


98 


1 


z 

H 
3 
O 
BS 
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C 


C                               «»  HAIN  PROGRAM  Ut  mi  1 

C  SAHI  2 

C  SAMI  3 

c     mimtmmmmmmnmtnmmtmmmmmtmmmsm  4 

C       n    SOLUTION  OF  SET  OF  N  4TH  ORDER  ORDINARY  DIFF,  EQUATIONS  »*SArtI  5 

C       »l  TO  DETERMINE  G(X)  FROM  WHICH  THE  STREAM  FUNCTION  ,  VELOCITY  tlSAMI  i 

C       n                            AND  PRESSURE  ARE  DEFINED  -««SAMI  7 

C       n                   (THIS  PROGRAM  WAS  WRITTEN  BY  L.  AMMAR)  **SAMI  8 

c     mmmtttmmtmmmmtmtmmmmtmmtxtmtmism  ? 

C  SAMI  10 

C          t   M  =  NUMBER  OF  POINTS  IN  THE  X-DIRECTION  SAMI  11 

C          *  N  =  NUMBER  OF  SET  OF  4TH  ORDER  EQUATIONS  (  INPUT  )  SAhi  12 

C           «   MN  =  DIMENSION  OF  GLOBAL  MATRIX  =  M»N  SAMI  13 

C          I   T  =  COEFFICIENTS  OF  MATRIX  OF  TYPE  1  SAHI  14 

C          »  R  =  COEFFICIENTS  OF  MATRIX  OF  TYPE  2  SAMI  15 

C          t  S'=  COEFFICIENTS  OF  MATRIX  OF  TYPE  3  SAMI  16 

C           »   RST  =  COEFFICIENTS  OF  GLOBAL  MATRIX  ^AMI  1? 

C          »  RR  =  INPUT  VECTOR  mi  is 

C           »   DX  =  STEP  SIZE  SAMI  19 

C          *  XAfXB  =  END  POINTS  ShNI  20 

C           »   X,HU  =  OBLATE  SPHEROIDAL  COORDINATES  SAM?  21 

C          t   ZETA  =  FIXES  THE  SIZE  OF  SPHEROIDAL  CAP  /  FEEDER  TUBE  SAHI  22 

C      •     t   RE  =  REYNOLDS  NUMBER  mi  23 

C          «  A.B»C»..,..Q  =  COEFFICIENTS  OF  DIFF,  EQUATIONS  SAMI  24 

C           *   WKAREA  =  WORK  AREA  OF  IMSL  OF  DIMENSION  ^  MN**2t3.*MN  SAHI  25 

C          *  PSI  =  STREAM  FUNCTION  SAMI  26 

C           *   W  =  DIMENSIONLESS  VELOCITY  (V»C«2/(NU*D))  SAMI  27 

C          I  P  =  DIMENSIONLESS  PRESSURE  (DEFINED  IN  THE  TEXT)  SAMI  23 

SAMI  29 

IMPLICIT  REALt8(A-H»0-Z)  mi  30 

DIMENSION  RST(303»303)fRSTINV(303,303),«KAREA(92718)rRR(303.M  SAMI  31 

A             6(303).GP(303)»GPP(303),GPPP(303)  qAMI 

COMMON/BLOKO/DX  SAi^I  33 

COHMON/BLOKl/A(101,3),B(101,3),C(101,3).D(101.31.E(10h3''  ^AMI  U 

C0MM0N/BL0K2/H(101..3),AI(101,3)fAJ(101,3)fAK(101.3)  mi 

COHMON/BLOK3./AL(10h3).AM(101.3),AN(101.3)..g(101,3)  mi  36 

COMMON/TYPET/TdOblOl)  SAMI  V 

COMMON/TYPER/R(101,101)  g^MI  ^3 

COMMON/TYPES/SdOi-lOl)  S^i^I  3. 

""^•^^  SAMI  40 
N=3 

MN=M*N 

XA=O.DO 


SAMI  41 

SAMI  42 

SAMI  43 

XB=1.D0  g^„j 

IRST=MN 
IDGT=6 
M1=M-1 


SAM:  45 
SAMI  46 
^•AMI  4"' 

DX=(XB-XA)/M1  3ft„j  ^8 


...INTRODUCE  THE  COEF.  OF  DIFF.  EQUATIONS...  SAHI  49 

DO  999  IPS=10,30,5  ^AHI  50 

DPS=IPS  SAhl  .1 

ZETA=DPS/10,  f "J  ^ 
DO  10  J=i!h 

DJ=J-1  ^^^^  5^ 

X=DJ»DX  ^^"^ 

Z=X»X+1./(ZETA»:ZETA)             .             *  SAHI  56 

"JtSr""  mil 

A(J»NI)=Z«Z»(X*X+(8.«iN»DN-4.»DN-3.)/((4.»DN+l.)»(4.ttN-3.)  SAMI  59 

B             »ZETA»ZETA))  ^^^^ 
B(J.NI)=-l6.»X<Z«iN*(2.»IiN-l.)/((4.*DNn.)»(4.ttiN-3.)«ZETA»ZETA)SAHI  61 

C(J»NI)=-4.»DN»(2.»DN-l.)*Z*(X«XH8.»IiN«iN-4.»DN-9.)  SAMI  62 

C             /U4.»DN+l.)»(4.*DH-3.)*ZETAtZETA))  SAHI  63 

D(JrNI)=8.»X»Z«iH*(2.*DH-l.)  SAMI  64 

E(J»NI)=4.HiN»(2.»DN-l.)»(DNn.)»(2.tIiN-3.)  SAMI  65 

D             »(Xm(B,»DN«iN-4.*IiN+l.)/((4.*DN+l.)»(4.»DN-3.)  SAMI  66 

E             tZETAIZETA))  SAMI  67 

H(J»NI)=Z»Z»2.»(DN-l.)*(2.«DN-3.)/((4.tIiN-5.)*(4,*DN-3.)  SAMI  68 

F             mZETAIZETA)  SAMI  69 

AI(J,NI)=8,mZ*(DN-l.)»(2.«DN-3.)/((4.«DN-5.)*(4.«DN-3.)  SAMI  70 

G               tZETA»ZETA)  SAMI  71 

AJ(JfNI)=-8.*Z»(DN-l.)»(DN-2.)»(2.«iN-3.)»(2,«DN-3.)  SAMI  72 

H               /((4.tDN-5.)»(4.«W-3,)«ZETA«ZETA)  SAMI  73 

AK(JfNI)=8.»(DN-l.)<(DN-l.)»(2.»DN-3.)»(2.*IiN-3,)»(IiN-2.)  SAMI  74 

I                 »(2.»DN-5,)/((4.IDN-5.)»(4,*DN-3.)»ZETA»ZETA)  SAMI  75 

AL( J»NI)=2.»Z»Z»(DN+1. )«(2.»DNn. )/( (4,»DK+1 . )*(4.*DN+3. )  SAMI  76 

J               »ZETA»ZETA)  SAMI  77 

AM(J>NI)=8.mZ»(DN+l.)»(2.«iNH.)/((4.*DN+l.)»(4.*DN+3.)  SAMI  78 

K               »ZETA»ZETA)  SAMI  79 

AN(J,NI)=-8.tZ»(DN+l.)»(DN+l.)»(2.»IiN+l.)*(2.»DN+3.)  SAMI  60 

L                     /((4,«DN+1.)»(4.»DN+3.)»ZETA»ZETA)  SAMI  81 

Q(J,NI)=8.»(DN+l.)*(DN+l.)»(DN+2.)t(2,»IiN+l.)*(2.*DN+l.)  SAMI  82 

M                »(2.tDN+3.)/((4,»DNn.)»(4,«DN+3.)»ZETA»ZETA)  SAMI  83 

CONTINUE  SAMI  84 

SAHI  85 

DO  999  IRE=lf3  SAMI  86 

DRE=IRE  SAHI  87 

RE=IiRE/10.  SAHI  88 

. . .CLEAR  GLOBAL  MATRIX  RST . . ,  SAMI  89 

SAHI  90 

DO  15  1=1 »MN  SAHI  91 

DO  15  .1=1  »HN  SAHI  92 

RST(IfJ)=0.D0  S**"!  93 

...COMPUTE  COEFFICIENTS  OF  MATRIX  T...  SAMI  94 

DO  30  I=1>N  SAHI  95 

CALL  TYPEKbM)  SAHI  96 

SAHI  97 

...ASSEHBLE  T  INTO  GLOBAL  HATRIX  RST...  SAHI  98 

DO  20  J=1»H  SAHI  99 

DO  20  K=1,M  SAMI 100 

JI=H»(I-1)+J  SAHIlOl 
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KI=M*(I-1)+K  SAHI102 

RST(JI.KI)=T(J>K)  SAMI103 

20    CONTINUE  SAHI104 

30    CONTINUE  SAMI105 

C  SAHI106 

C       ...COMPUTE  COEFFICIENTS  OF  rtATRIX  R...  SAMI107 

DO  50  I=2>N  SAMI108 

CALL  TYP£2(I»M)  SAHI109 

C  SAHIllO 

C       ...ASSEMBLE  R  INTO  GLOBAL  MATRIX  RST...  SAMIIU 

DO  40  J=1jM  SAhI112 

DO  40  K=1tM  SAMI113 

JI=H»(I-1)+J  SAMI114 

KI=Mt(I-2)+K  SAMI115 

RST(JIfKI)=R(J»K)  SAMI116 

40     CONTINUE  SAMI117 

50    CONTINUE  SAMI118 

C  SAMIll? 

C       ...COMPUTE  COEFFICIENTS  OF  MATRIX  S...  SAMI120 

N1=N-1  SAMI121 

DO  70  I=1»N1  SAMI122 

CALL  TYPE3(IfM)  SAMI123 

C  SAMI124 

C       ...ASSEMBLE  S  INTO  GLOBAL  MATRIX  RST,, .  SAMI125 

DO  60  J=1..M  SAMI126 

DO  60  K=1»M  SAMI127 

JI=Mt(I-l)+J  SAMI128 

KI=M*I+K  3AMI129 

RST(JIfKI)=S(JfK)  SAMI130 

60    CONTINUE  SAMI131 

70    CONTINUE  SAMI132 

C  SAMI133 

C       ...CLEAR  INPUT  VECTOR  RR,,,  SAHI134 

DO  80  I=hMN  SAMI135 

80        RR(I)=O.DO  SAMI136 

C       ...INTRODUCE  DIRICHLET  BC  AT  XA...  SAMI137 

MM«=(N-1)*M+1  SAMI138 

DO  94  J=hMMM.H  SAMI139 

G(J)=0,DO  SAMI140 

G(l)=3./10.tRE  SAMI141 

G(M+1)=1./30.»RE  SAMI142 

IF(G(J).NE,0,)  GO  TO  91  SAMI145 

DO  90  K=1,MN  SAMI144 

RST(JfK)=0.DO  SAMI145 

90  RST(K,J)=0.D0  SAMI146 
RST(JjJ)=l,DO  3AHI147 

GO  TO  94  SAMI148 

C  SAMI149 

C       ...COMPUTE  INPUT  VECTOR  RR(I)  MODIFIED  BY  BC  AT  XA..,  SAMI150 

91  CONTINUE  SAMI151 
DO  92  KX=1,HN  SAMil52 

92  RR(KX)=RR(KX)-RST(KX!J)»G(J)  SAhI153 
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DO  93  KX=hMN  3AI1I154 

RST(KXfJ)=O.DO  SAiiI155 

RST(J,KX)=O.DO  SArtI156 

93  CONTINUE  SAIiI157 

RST(J»J)=l.IiO  SAHI158 

RR(J)=G(J)  SAMIIS'? 

94  CONTINUE  SAMI160 
C  SAI1I161 
C       ...INTRODUCE  DIRICHLET  BC  AT  XB...  SAMI162 

DO  96  -MmrH  smibz 

G(J)=O.DO  SAHI164 

DO  95  K=1»HN  SAHI165 

RST(J>K)=O.DO  SAMI166 

95  RST(K»J)=0.DO  SAHI167 
RST(J>J)=1,D0  SAMI168 

96  CONTINUE  3AMI169 
C  SAHI170 
C  ...  INVERT  GLOBAL  HATRIX  RST  AND  hULTIPLY  BY  VECTOR  RR...  3AliI171 
C       ...USING  IMSL  LIBRARY.,.  SAMI172 

CALL  LINV2F(RSTfHN.IRST.RSTINV.IDGTf«KAREATlER)  SAHI)73 

DO  97  I=1»HN  SAhI174 

G(I)=O.DO  3AiiI175 

DO  97  J=1jHN  SAhI176 

97  6(I)=6(I)+RSTINV(I.J)»RR(J)  3AHI177 
C  SAHI173 
C       ...COMPUTE  DERIVATIVES  OF  G...  3AMI179 

DO  98  J=hNHH.N  SAhllbO 

GP(J)=O.DO  SAHIiai 

GPP(J)=2./(IiX»DX)l(G(J+l)-G(J))  SAt1IlS2 

GPPP(J)=l./(2.*DXtDXtDX)t(-5.*G(J)+13.«G(Jtl)-24.»G(J+2)  SAHI1S-. 

N               H4,»6(J+3)-3.»G(Jt4))  3AMI134 

98  CONTINUE  3AHI155 
DO  99  J=«jHNfH  3AHI186 

GP(J)=O.DO  SAhI187 

GPP(J)=2,/(DXtDX)t(G(J-l)-G(J))  SAHI188 

GPPP(J)=l./!2.»DX*DX*DX)*(5.*G(J)-i8.*G(J-l)+24,*6(J-2)  SAHI189 

Q               -14.»G(J-3)+3.»G(J-4))  SArillyO 

99  CONTINUE  SAMIl?! 
DO  101  J=2»i11  SAnll?: 
«J=(N-i;itH+J  SAM  1 193 
DO  101  JX=J»HJ)H  SHrtIi"T4 

6P(JX)=1./(2.*DX)»(G(JX+1)-GUX-1))  SAnll?5 

101  GPP(JX)=l./(DX«iX)*(GUX-l/-2.»G(JX)+G(JX+l))  3AhIlV6 
LJ=(N-l)tH+2  SAfillv? 
DO  102  JX=2»LJ»H  SftCll98 

GPPP  ( JX ) = 1 ,  /  ( 2 .  *DX*DX»DX )  t  ( -G  ( JX )  +2 .  *6  ( JX  - 1 )  -2 .  *G  ( JX+ 1 )  3  Ali  11  ?  9 

R                +G(JX+2))  SAhlZOO 

102  CONTINUE  3Afii:0i 
HNl=MN-i  SAHi::2 
DO  103  JX=«l,HNliH  SAhI203 

GPPP(JX)=l./(2.*DX«DX«DX)l(GiJX)-2.«6(JX+l)f2.»G(JX-l)-6<JX-2))5A«I204 

103  CONTINUE  SA«I205 
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M2=H-2 

SAHI206 

DO  104  J=3fH2 

SAMI207 

HJ=(N-l)»h+J 

SAMI208 

DO  104  JX=J»hJ»M 

SAMI209 

GPFf(JX)=l,/(2,*DX*DXMX)«<-6(JX-2)+2.»G(JX-l)-2.*G(JX+l) 

SAMI210 

i  +G(JX+2)) 

SAMI211 

CONTINUE 

SAMI212 

HRITE(6»222)NiZETA»RE 

SAMI213 

DO  111  IPR=1,H 

SAhI214 

DJ=IPR-1 

SAMI215 

X=DJ*DX 

SAMI216 

URITE(6(333)X 

SAMI217 

IIP=(N-l)tH+IPR 

SAMI218 

KFR=-1 

SAMI219 

T\n    <  4  <      ir>f\    Trir\    ttpi  vi 

DO  111  JPR=IFR>IIP»M 

SAMI220 

KPR=KPR+2 

SAMI221 

WRITE(i.444)KPR»G(JPR),KPRj6P(JPR)»KPR»GPP(JPR),KFR»6PPP(JPR) 

SAMI222 

CONTINUE 

SAMI223 

«RITE(6»555)IER 

SAMI224 

SAMI225 

...COMPUTATION  OF  THE  STREAM  FUNCTION  (DIMENSIONLESS) . . . 

SAMI226 

URITE(6,888) 

SAMr227 

DO  105  J=1»M 

SAMI228 

DJ=J-1 

SAMI229 

X=DJ»DX 

SAMI230 

WRITE(6>666)X 

SAHI231 

DO  105  K=1j11 

SAMI232 

DK=K-1 

SAMI233 

SAMI234 

1 II 1— 1 141 1 

uu=u»u 

SAMI235 

UUUU=UU»UU 

SAMI236 

PSI=(1.-UU)»(G(J)+G(H+J)«(5.*UU-1,)«,/2.IG(2»H+J) 

SAMI237 

t(21.*UUUU-14,«UU+l,)»15,/8.) 

SAMI233 

URITE(6»777)U>PSI 

SAMI239 

SAMI240 

...COMPUTATION  OF  SPEED  DISTRIBUTION... 

SAMI241 

URITE(6»889) 

3AMI242 

DO  107  J=1»H 

SAHI243 

DJ=J-1 

3AMI244 

X=DJ*DX 

SAMHtS 

yRITE(6j666)X 

SHfiI2'>6 

DO  107  K=hll 

SAhI247 

DK=K-1 

SAMI248 

U=DK/10. 

SAi1I249 

UU=U*U 

SAMI250 

UUU=UUtU 

SAMI251 

UUUU=UUU»U 

3ArtI252 

UUUUU=UUUU»U 

3AMI253 

Z=XmUU/(ZETA»ZETA) 

SAMI254 

ZZ=X*X+1./(ZETA»ZETA) 

ShMI255 

IF(X.NE.0.OR.U.NE.O.;  30  TO  106 

SAMI256 

H=O.DO 

SAhI257 

GO  TO  107 

SHr1I25S 
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106  W=l , / DSQRT ( Z ) »DSQRT  ( 4 . /lit ( G ( J ) »U+3 . *G ( H+ J ) * ( 3 . WUU-3 , *U )  SAMI259 
U  +15./8.tG<2»H+J)»(63,»UUUUU-70,*UUU+15,*U))«2  SAMI260 

V  +(l.-UU>»(GP(J)+3./2.*6F(«IJ)»(5,*UU-l.)  3AhI261 
U  +15,/8.«GP(2»M+J)»<21.«UUUU-14.»UU+l.))»t2)  SAHI262 

107  WRITE(6r890)U.«  SAHI263 
C  3AMI2A4 
C       ...COMPUTATION  OF  THE  PRESSURE.,.  SAMI265 

wRiTE(6»89i)  mim 

DO  108  J=2»H  SAI1I2d7 

DJ=J-1  SAHI26e 

x=DJ*[ix  mim 

URITE(6.666)X  3AMI270 

ZZ=Xmi./(ZETA»ZETA)  3Ai1I271 

61=ZZ«GPP(J)-2.tG(J)  3AhI272 

GP 1 =ZZ»GPPP ( J ) +2 . »X*GPP ( J ) -2 . *GP ( J )  3 AM I 273 

G3=ZZ»GPP(M+J)-12.»G(M+J)  SAMI274 

6P3=ZZ»GPPP  ( M+ J )  +2 .  »X»GPP  ( H+ J )  - 12 .  »GP  ( M+ J )  SAMI275 

G5=ZZ*GPP(2*M+J)-30.»G(2*M+J)  SAMI276 

GP5=ZZ»GPPP  ( 2tM+ J)  +2 .  »X»GPP  ( 2»M+ J )  -30 .  tGP  ( 2*M+ J )  3  AMI277 

Zl =1 . / ( X*X ) -5 . »ZETA*ZET A  SAMI278 

Z2=63 .  *X«X»ZETA*ZETA»ZETA*ZETA-14 . »ZETA»ZETA-1 . / ( X»X )  SAMI279 

Z3=l./X+5.»XtZETA»ZETA  SAMI280 

Z4=21.*X»XtX»ZETAtZETA*ZETA«ZETA+14.*ZETAXZETA»X+l,/X  3AHI281 

DO  108  K=l»ll  SAMr232 

DK=K-1  SAMI283 

U=riK/10.  SAMI2S4 

UU=U*U  3AMI285 

Z=XIX+UU/(ZETA#ZETA)  SAMI286 

P=(-Gl/(X»X)+3./2.»Zl*G3+15./8.tZ2*G5+GPl/X-3./2.*Z3»6P3  SAMI287 

X   +15./8.«4»GP5)»ZETA«DATAN(U/(X»ZETA))  SAMI288 

Y  +(-Gl/X+3./2.IZ3»G3+15./8.t(21.IX*UUtZETA«ZETA.14.mZETA*ZETA  3AMI239 
Z  -l./X)»G5)tU/Z+(GP3+(21.l(UU/3.-X»X*ZETA«ZETA)-14,)*GP5/4.  SAMI290 
A   -63./4.»X»ZETAtZETA*G5)tl5,/2,»U*ZETA»ZETA  3AMI291 

108  WRITE(6.892)U>P  3AMI292 
999   CONTINUE  3AMI293 

STOP  SAM I 294 

222   FORMAT  (1  HI , 5X . '  N= ' ..  12  j  5X » ' ZETA= ' » F5 . 3 1 5X » '  RE= ' . F5 . 3 ..  /// )  S AHI 295 

333  F0RMAT(3X»'X='»F4.2)  z^nVlih 
444   F0RMATa4X.'6M2j'  =  '.E11.4.5X»'6P'»I2»'  =  'iE11.4.5XrGpp',I2..  ='.ESAMI297 

Ull,4»5Xj'GPPP'fI2»'=')E11.4j/)  5ArtI298 

555   F0RMAT(5X»'IER='»I3.////)  SAMI299 

666   F0RMAT(3Xf'X='rF4.2)  SAMI300 

777   F0RMAT(20X)'rtU='.F4.2.5X>'PSI='»E15.5T./)  SAMI301 

888  F0RHAT(30X)'DIMENSI0NLES3   STREAM  FUNCTION   FSKX-MU)  ',///)  SAMI302 

889  F0RMAT(lHl»30Xf 'SPEED  DISTRIBUTION  y*Ctt2/(NU»D)'i///)  SAMI303 

890  F0RMAT(20Xf'MU='fF4.2j5X>'y='.E15,5f./)  SAhI304 

891  F0RMAT(lHli50X» 'PRESSURE  DISTRIBUTION'.///)  SAMI305 

892  F0RMAT(20X»'MU=  ,F4.2.'P='iE15.5j/)  •3AMI306 

END                                                           •  3AMI307 
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Q  SAMI303 

C                                   m  TYPEl  nt  SArtI309 

C  SAHI310 

SUBROUTINE  TYPEl (L»M)  3A«I311 

IMPLICIT  REALt8(A-H»0-Z)  SAh!512 

COHMON/BLOKO/DX  3AhI313 

CO«HON/BLOKl/A(101»3)»B(101j3))C(10b3)fD(101.3)»E(10b3)  SAHI3i4 

COMMON/TYPET/TdOblOl)  3AHI315 

C       .,, CLEAR  MATRIX  T...  SAMI3U 

DO  100  1=1. M  3AMI317 

DO  100  J=1>H  SAMI318 

100       T(IfJ)=O.DO  SAMI319 

C       . . .COMPUTE  T  FOR  NI=L . . ,  3AMI320 

T(2»2)=7.»A(25L)-B(2.L)»DX/2.-2,*C(2»L)*DX*DX  SAhI321 

1  +E(2>L)tDX»DXtDX*DX  3AHI322 
T(M-hM-l)=7.*A(M-hL)+B(M-l.L)*DX/2.-2.»C(H-lfL)IDX«iX  3AhI323 

2  +E(M-bL)«DX»DXtDX«iX  SAMI324 
M2=M-2  SHflI325 
DO  110  I=3»M2  SAMI326 

110       T ( I , I ) =6 . t A( I » L)-2 , *C ( I J  L ) *DX»DX+E (I » L ) »DX*DX*DX»DX  SAMI327 

H1=H-1  3AHI32S 

DO  120  I=2»M1  SAhI329 
T(  b I -1 ) =-4 . tA (I » L ) +B  < I >  L ) »DX+C (I >  L ) *DX»DX-D ( 1 1 L ) *DX*DX*DX/2 .  SAMI330 

120       T(hI+l)=-4.»A(I»L)-B(I.L)»DX+C(I.L)*DX»DXKi(I»L)»DX-.HiX*DX/2.  3A(iI331 

DO  130  I=2jH2  SArtI332 

T(HljI-l)=A(I+l»L)-B(in.L)*DX/2.  ShMI333 

130       T(I»I+2)=A(I.L)+B(IjL)»DX/2.  SArtI334 

RETURN  SAM I 335 

END  SAMI336 
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C  SAMI337 

C                                   tn  TYPE2  m  SAHI338 

C  SAHI33? 

SUBROUTINE  TYPE2(L»H)  SAHI340 

IMPLICIT  REAL»8(A-H»0-Z)  SAMI341 

COM«ON/BLOKO/DX  SAi1I342 

C0HM0N/BL0K2/H(101»3)»AI(101j3)»AJ(10h3)»AK(101i3)  SA«I343 

COHHON/TYPER/R(101»101)  SAi1I344 

C       ...CLEAR  HATRIX  R...  SAMI345 

DO  200  1=1, M  3AHI346 

DO  200  J=1,M  SAHI347 

200       R(I,J)=0,D0  SAMI348 

C       ...COhPUTE  R  FOR  NI=L...  SAHI349 

R(2,2)=7.*H(2,L)-flI(2»L)»DX/2.-2.»AJ(2,L)»DX»DX  SAMI350 

3  +AK(2fL)»DXtDX»DX*DX  SAMI351 
R(M-l»H-l)=7.«H(H-l,L)+AI(H-l,L)»DX/2.-2.»AJ(«-l.L)*DX«DX  SAMI352 

4  +AK(M-l,L)*DX«DXtDX*DX  SAHI353 
M2=M-2  3AMI354 
DO  210  1=3, H2  SAMI355 

210       R(I,I)=6.tH(I,L)-2,»AJ(I>L)tDXtDX+AK(I,L)*DX*DXtDX«iX  SAI1I356 

M1=H-1  SAMI357 

00  220  1=2, HI  SriHI358 

R(I,I-l)=-4.»H(I,L)+AI(I,L)»X+AJ(I,L)tDX»DX  3AhI359 

220       R(I,I+1)=-4.»H(I,L)-AI(I,L)»DX+AJ(I,L)»DX»DX  SAiiI360 

DO  230  1=2, M2  SAHI361 

R(I+l,I-l)=H(I+l,L)-AI(I+l,L)»DX/2.  SAMI3i2 

230       R(I,I+2)=H(I»L)+AI(I,L)»DX/2,  3AHI3o3 

RETURN                                                          ■          '  .  SAHI3i4 

END                                                                     ■•  SAhI365 
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C  SAMI366 

C                                   m  TYPE3  m  SArtI367 

C  3AflI368 

SUBROUTINE  TYPE3(L,M)  SAHI369 

IHPLICIT  REAL»8(A-H»0-Z)  SAIiI370 

COMhON/BLOKO/DX  SAHI371 

C0HH0N/BL0K3/AL(  101  »3) » AM(  101, 3) » AN(  101 ,3) » Q(  10b3)  SAhI372 

COHHON/TYPES/SdOhlOl)  SA(1I373 

C       ...CLEAR  HATRIX  S...  SAMI374 

DO  300  1=1, H  SAMI375 

DO  300  J=hH  SAMI376 

300       S(I,J)=O.DO  3AMI377 

C       ...COMPUTE  S  FOR  NI=L...  SAMI378 

S(2,2)=7.*AL(2,L)-AH(2,L)»DX/2.-2.»AN(2,L)tDXtDX  SAMI379 

5  +Q(2»L)HiX«iXiDX«iX  SAHI380 
S(M-i,M-l)=7.«AL(M-l,L)+AM(M-l,L)*DX/2.-2.*AN(M-l,L)*DX«iX  SAMI331 

6  +Q(M-l,L)»nX»DX»DX»DX  SAHI382 
M2=H-2  SAMI383 
DO  310  I=3.H2  SAhI384 

310       S ( I , I ) =6 . lAL { I, L ) -2 , »AN (I , L ) »DXIDX+Q (I , L ) tDX»DX«DX»DX  SAMI385 

'<1=«-1  SAHI386 

DO  320  I=2,H1  SAHI387 

S(I,I-1)=-4.IAL(I,L)+AH(I,L)«DX+AN<I,L)*DX»DX  SAHI388 

320       S(I,I+l)=-4.*AL(I,L)-AH(I,L)!»DX+AN(I,L)«iXIDX  3AMI389 

DO  330  I=2,H2  3AHI390 

S(IH,I-l)=AL(Itl,L)-AM(I+l,L)»DX/2.  SAHI391 

330       S(I,I+2)=AL(I,L)+A«(I,L)«DX/2.  SAMI392 

RETURN  smm 


SAMI394 


108 


A  Sample  of  the  Computer  Output 


ZETA=2,000  RE=0.100 


G  1=  0.3000D-01 
G  3=  0,3333D-02 
G  5=  0.0 

6  1=  0,30000-01 
6  3=  0.3330D-02 
G  5=-0.7708D-09 

G  1=  0,2998D-01 
G  3=  0.3321D-02 
G  5=-0.2826D-08 

G  1=  0.2996D-01 
G  3=  0.3307D-02 
G  5=-0.5956D-08 

G  1=  0.2993D-01 
G  3=  0,3287D-02 
G  5=-0.9997D-08 

G  1=  0.2990D-01 
G  3=  0.3264D-02 
G  5=-0.1481D-07 

G  1=  0.2986D-01 
6  3=  0,32360-02 
G  5=-0. 20300-07 

G  1=  0,29810-01 
G  3=  0,32060-02 
G  5=-0, 26330-07 

G  1=  0,29760-01 
G  3=  0,31720-02 
0  5=-0. 32990-07 

G  1=  0.29700-01 
G  3=  0.31350-02 
G  5=-0. 40070-07 

G  1=  0.29640-01 
G  3=  0.30960-02 
G  5=-0. 47610-07 


GP  1=  0,0 
GP  3=  0,0 
GP  5=  0,0 

GP  l=-0, 91880-03 
GP  3=-0, 61 740-03 
GP  5=-0, 14130-06 

GP  l=-0, 17430-02 
GP  3=-0. 11770-02 
6P  5=-0, 25930-06 

GP  l=-0,2488D-02 
GP  3=-0, 16820-02 
GP  5=-0, 35860-06 

GP  l=-0, 31660-02 
GP  3=-0, 21370-02 
GP  5=-0. 44280-06 

GP  l=-0. 37900-02 
GP  3=-0, 25460-02 
GP  5=-0, 51530-06 

GP  l=-0,43730-02 
GP  3=-0. 29120-02 
GP  5=-0, 57840-06 

GP  l=-0, 49240-02 
GP  3=-0, 32380-02 
GP  5=-0. 63430-06 

GP  l=-0. 54540-02 
GP  3=-0. 35280-02 
GP  5=-0, 68470-06 

GP  l=-0, 59720-02 
6P  3=-0, 37340-02 
GP  5=-0. 73090-06 

GP  l=-0, 64840-02 
GP  3=-0. 40 100-02 
GP  5=-0. 77400-06 


GPP  l=-0.96990-01 
GPP  3=-0,6475D-01 
GPP  5=-0. 15420-04 

GPP  l=-0, 86780-01 
GPP  3=-0. 58730-01 
GPP  5=-0, 12840-04 

GPP  l=-0, 78080-01 
GPP  3=-0, 53130-01 
GPP  5=-0, 10760-04 

GPP  l=-0, 70800-01 
GPP  3=-0, 47930-01 
GPP  5=-0, 90930-05 

GPP  l=-0. 64830-01 
GPP  3=-0, 431 10-01 
GPP  5=-0. 77660-05 

GPP  l=-0, 60080-01 
GPP  3=-0, 38650-01 
GPP  5=-0, 67200-05 

GPP  l=-0, 56450-01 
GPP  3=-0. 34530-01 
GPP  5=-0. 59040-05 

GPP  l=-0, 53840-01 
GPP  3=-0, 30730-01 
GPP  5=-0, 52750-05 

GPP  l=-0. 52160-01 
GPP  3=-0, 27230-01 
GPP  5=-0, 47980-05 

GPP  l=-0, 51310-01 
GPP  3=-0, 24020-01 
GPP  5=-0. 44430-05 

GPP  l=-0. 51220-01 
GPP  3=-0. 21070-01 
GPP  5=-0. 41860-05 


GPPP  1=  0,10820+01 
GPPP  3=  0,61960+00 
GPPP  5=  0,26980-03 

GPPP  1=  0,94560+00 
GPPP  3=  0,58080+00 
GPPP  5=  0,23280-03 

GPPP  1=  0,79890+00 
GPPP  3=  0,53990+00 
GPPP  5=  0.18740-03 

GPPP  1=  0.66230+00 
GPPP  3=  0.50110+00 
GPPP  5=  0,14970-03 

GPPP  1=  0.53580+00 
GPPP  3=  0,46410+00 
GPPP  5=  0,11860-03 

GPPP  1=  0,41910+00 
GPPP  3=  0.42920+00 
GPPP  5=  0.93100-04 

GPPP  1=  0.31220+00 
GPPP  3=  0,39610+00 
GPPP  5=  0,72250-04 

GPPP  1=  0,21460+00 
GPPP  3=  0,36490+00 
GPPP  5=  0.55310-04 

GPPP  1=  0,12630+00 
GPPP  3=  0,33550+00 
GPPP  5=  0,41610-04 

GPPP  1=  0,46860-01 
GPPP  3=  0.30790+00 
GPPP  5=  0.30600-04 

GPPP  l=-0. 24010-01 
GPPP  3=  0.28200+00 
GPPP  5=  0.21790-04 
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X=0.11 


X=0.12 


X=0.13 


X=0,14 


X=0,15 


X=0.16 


X=0.17 


X=0.18 


X=0.19 


X=0,20 


X=0.21 


x=o.: 


X=0.23 


G  1=  0.2957D-01 
G  3=  0.3055D-02 
G  5=-0.5555D-07 

6  1=  0.2950D-01 
G  3=  0.3012D-02 
G  5=-0.6390D-07 

G  1=  0,2942D-01 
G  3=  0.2967D-02 
G  5=-0.72A4D-07 

G  1=  0,29340-01 
G  3=  0.2921D-02 
G  5=-0,8177D-07 

G  1=  0.2925D-01 
G  3=  0.2874D-02 
G  5=-0,9127D-07 

6  1=  0.2915D-01 
G  3=  0.2826D-02 
G  5=-0.1011D-06 

G  1=  0.2905D-01 
G  3=  0.2777D-02 
G  5=-0.1114D-06 

G  1=  0.2894D-01 
G  3=  0.2728D-02 
G  5=-0.1220D-06 

G  1=  0. 288311-01 
G  3=  0.2678D-02 
G  5=-0,1331D-06 

G  1=  0.2871D-01 
G  3=  0,26270-02 
G  5=-0.1445D-06 

G  1=  0.2858D-01 
G  3=  0,2576D-02 
G  5=-0,1563D-06 

G  1=  0.2844D-01 
G  3=  0.2525D-02 
G  5=-0.1685D-06 

G  1=  0.2830D-01 
G  3=  0.2474D-02 
G  5=-0.1B12D-06 


GP  l=-0.6999D-02 
GP  3=-0,4207D-02 
GP  5=-0.8150D-06 

GP  l=-0,7523D-02 
GP  3=-0,4378D-02 
GP  5=-0,8545D-06 

GP  l=-0,8061D-02 
GP  3=-0,4526D-02 
GP  5=-0.8930D-0A 

GP  l=-0,8618D-02 
GP  3=-0.4653D-02 
GP  5=-0.9311D-06 

GP  l=-0.9197D-02 
GP  3=-0.4760D-02 
GP  5=-0,9689D-06 

GP  l=-0,9803D-02 
GP  3=-0.4849D-02 
GP  5=-0,1007D-05 

GP  1=-0,1044D-01 
GP  3=-0,4923D-02 
GP  5=-0.1045D-05 

GP  1=-0.1110D-01 
GP  3=-0,4982D-02 
GP  5=-0.1084D-05 

GP  1=-0,1180D-01 
GP  3=-0,5028D-02 
GP  5=-0,1123D-05 

GP  1=-0,1254D-01 
GP  3=-0.5062D-02 
GP  5=-0.1162D-05 

GP  1=-0.1330D-01 
GP  3=-0,5086D-02 
GP  5=-0,1202D-05 

GP  1=-0,1410D-01 
GP  3=-0.5101D-02 
GP  5=-0.1242D-05 

GP  1=-0.1494D-01 
GP  3=-0.5107D-02 
GP  5=-0. 128311-05 


GPP  1=-0,5179D-01 
GPP  3=-0.1838D-01 
GPP  5=-0,4007D-05 

GPP  1=-0.5295D-01 
GPP  3=-0.1592D-01 
GPP  5=-0,3890D-05 

GPP  1=-0,5462D-01 
GPP  3=-0.1367D-01 
GPP  5=-0,3821D-05 

GPP  l=-0. 56730-01 
GPP  3=-0.1163D-01 
GPP  5=-0.3789D-05 

GPP  l=-0, 59200-01 
GPP  3=-0,9783D-02 
GPP  5=-0, 37860-05 

GPP  l=-0, 61980-01 
GPP  3=-0. 81060-02 
GPP  5=-0, 38030-05 

GPP  1 =-0,65000-01 
GPP  3=-0. 65910-02 
GPP  5=-0. 38350-05 

GPP  l=-0, 68210-01 
GPP  3=-0, 52240-02 
GPP  5=-0. 38760-05 

GPP  1=-0.7155D-01 
GPP  3=-0. 39950-02 
GPP  5=-0.39220-05 

GPP  l=-0, 74970-01 
GPP  3=-0. 2891 0-02 
GPP  5=-0, 39700-05 

GPP  l=-0, 78440-01 
GPP  3=-0. 19030-02 
GPP  5=-0. 40140-05 

GPP  l=-0, 81900-01 
GPP  3=-0, 10210-02 
GPP  5=-0. 40540-05 

GPP  l=-0. 85320-01 
GPP  3=-0. 23460-03 
GPP  5=-0. 40860-05 


GPPP  l=-0. 86670-01 
GPPP  3=  0.25780+00 
GPPP  5=  0.14810-04 

GPPP  l=-0. 14150+00 
GPPP  3=  0.23520+00 
GPPP  5=  0.9313D-05 

GPPP  l=-0. 18880+00 
GPPP  3=  0,21410+00 
GPPP  5=  0.50340-05 

GPPP  l=-0. 22900+00 
GPPP  3=  0.1945D+00 
GPPP  5=  0.17560-05 

GPPP  l=-0. 26260+00 
GPPP  3=  0,17640+00 
GPPP  5=-0.6952D-06 

GPPP  l=-0,2899D+00 
GPPP  3=  0,15960+00 
GPPP  5=-0. 24620-05 

GPPP  l=-0, 31 130+00 
GPPP  3=  0,14410+00 
GPPP  5=-0. 36580-05 

GPPP  l=-0. 32730+00 
GPPP  3=  0.12980+00 
GPPP  5=-0. 43700-05 

GPPP  l=-0. 33820+00 
GPPP  3=  0.11660+00 
GPPP  5=-0. 46670-05 

GPPP  1 =-0.34440+00 
GPPP  3=  0.10460+00 
GPPP  5=-0. 46040-05 

GPPP  1 =-0,34630+00 
GPPP  3=  0.93530-01 
GPPP  5=-0. 42230-05 

GPPP  l=-0.3443D+00 
GPPP  3=  0.83440-01 
GPPP  5=-0. 35560-05 

GPPP  l=-0. 33860+00 
GPPP  3=  0.74230-01 
GPPP  5=-0, 26290-05 


G  1=  0.2814D-01 
G  3=  0.2423D-02 
G  5=-0,1942D-06 

G  1=  0.2798D-01 
G  3=  0,2372D-02 
G  5=-0,2076D-06 

G  1=  0,2781D-01 
G  3=  0.2321D-02 
G  5=-0,2215D-06 

G  1=  0,2763D-01 
G  3=  0.2270D-02 
G  5=-0.235BD-06 

G  1=  0.2744D-01 
G  3=  0,2220D-02 
G  5=-0.2504D-06 

G  1=  0.2724D-01 
6  3=  0,2170D-02 
G  5=-0,2655D-06 

G  1=  0,2702D-01 
G  3=  0.2120D-02 
G  5=-0.2810D-06 

G  1=  0,2680D-01 
G  3=  0.2070D-02 
G  5=-0.2969D-06 

G  1=  0,2(S57D-01 
G  3=  0.2021D-02 
G  5=-0.3131D-06 

G  1=  0.2633D-01 
G  3=  0.1972D-02 
G  5=-0,3297D-06 

G  1=  0,2407D-01 
6  3=  0,1923D-02 
G  5=-0.3467D-06 

G  1=  0.2580D-01 
G  3=  0.1875D-02 
6  5=-0,3640D-06 

G  1=  0,2553D-01 
G  3=  0,1827D-02 
G  5=-0.3816D-06 


GP  1=-0.1581D-01 
GP  3=-0.5106D-02 
GP  5=-0.1324D-05 

GP  1=-0.1671D-01 
GP  3=-O.5098D-O2 
GP  5=-0.1365D-05 

GP  1=-0.1765D-01 
GP  3=-0.5085D-02 
GP  5=-0.1406D-05 

GP  1=-0.1861D-01 
GP  3=-0.5066D-02 
GP  5=-0.1447D-05 

GP  1=-0.1961D-01 
GP  3=-0,5043D-02 
GP  5=-0.1488D-05 

GP  1=-0.2063D-01 
GP  3=-0.5016D-02 
GP  5=-0.1528D-05 

GP  1=-0.2167D-01 
GP  3=-0.4985D-02 
GP  5=-0.15&7D-05 

GP  1=-0,2274D-01 
GP  3=-0.4951D-02 
GP  5=-0.1606D-05 

GP  1=-0.2383D-01 
GP  3=-0,4915D-O2 
GP  5=-0.1643D-05 

GP  1=-0,2493D-01 
GP  3=-0,4877D-02 
GP  5=-0.1678D-05 

GP  1=-0.2605D-01 
GP  3=-0,4836D-02 
GP  5=-0.1712D-05 

GP  1=-0,2719D-01 
GP  3=-0.4794D-02 
GP  5=-0.1744D-05 

GP  1=-0.2833D-01 
GP  3=-0.4751D-02 
GP  5=-0.1773D-05 


GPP  1=-0.88671)-01 
GPP  3=  0.4638D-03 
GPP  5=-0,4107D-05 

GPP  1=-0.9192D-01 
GPP  3=  0.1083D-02 
GPP  5=-0,4115D-05 

GPP  1=-0.9503D-01 
GPP  3=  0,1629D-02 
GPP  5=-0.4108D-05 

GPP  1=-0.9799D-01 
GPP  3=  0,2110D-02 
GPP  5=-0.4084D-05 

GPP  1=-0.1008D+00 
GPP  3=  0,2533D-02 
GPP  5=-0.4042D-05 

GPP  l=-0.1033D+00 
GPP  3=  0.2902D-02 
GPP  5=-0.3978D-05 

GPP  l=-0.1057D+00 
GPP  3=  0.3224D-02 
GPP  5=-0.3892D-05 

GPP  1=-0.1078D+00 
GPP  3=  0.3503D-02 
GPP  5=-0.3782D-05 

GPP  1=-0,1097D+00 
GPP  3=  0.3745D-02 
GPP  5=-0.3645D-05 

GPP  1=-0,1114D+00 
GPP  3=  0,3953D-02 
GPP  5=-0.3482D-05 

GPP  l=-0.1127D+00 
GPP  3=  0.4131D-02 
GPP  5=-0.3290D-05 

GPP  l=-0.1138D+00 
GPP  3=  0.4283D-02 
GPP  5=-0.3068D-05 

GPP  1=-0.1146D+00 
GPP  3=  0.4413D-02 
GPP  5=-0.2815D-05 


GPPP  l=-0.3297D+00 
GPPP  3=  0,6584D-01 
GPPP  5=-0,1462IH)5 

GPPP  1=-0.3179D+00 
GPPP  3=  0.5826D-O1 
GPPP  5=-0.7167D-07 

GPPP  1=-0.3035D-H)0 
GPPP  3=  0.5139D-01 
GPPP  5=  0.1529D-05 

GPPP  1=-0.2867IH00 
GPPP  3=  0.4518D-01 
GPPP  5=  0.3328D-05 

GPPP  l=-0.2678D+00 
GPPP  3=  0.3958D-01 
GPPP  5=  0,5315D-05 

GPPP  1=-0,2470D+00 
GPPP  3=  0.3456D-01 
GPPP  5=  0.7482D-05 

GPPP  l=-0,2247D+00 
GPPP  3=  0.3006D-01 
GPPP  5=  0.9820D-05 

GPPP  l=-0.2009D+00 
GPPP  3=  0,2605D-01 
GPPP  5=  0.1232D-04 

6PPP  l=-0,1760D+00 
GPPP  3=  0.2248D-01 
GPPP  5=  0.1497D-04 

GPPP  l=-0.1501D+00 
GPPP  3=  0.1931D-01 
GPPP  5=  0.1777D-04 

GPPP  l=-0.1233D+00 
GPPP  3=  0.1652D-01 
GPPP  5=  0.2070D-04 

GPPP  1=-0.9575D-01 
GPPP  3=  0,1407D-01 
GPPP  5=  0.2375DH)4 

GPPP  1=-0.6768D-01 
GPPP  3=  0.1193D-01 
GPPP  5=  0.2691D-04 
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X=0,37 


X=0,38 


X=0.39 


X=0,40 


X=0,41 


X=0.42 


X=0,43 


X=0.44 


X=0,45 


X=0.46 


X=0,47 


X=0,48 


X=0.49 


G  1=  0,2524D-01 
G  3=  0.1730D-02 
G  5=-0.3994D-06 

6  1=  0.2494D-01 
G  3=  0.1733D-02 
G  5=-0.4176D-06 

G  1=  0.24A3D-01 
G  3=  0.1687D-02 
G  5=-0, 43590-06 

G  1=  0.2430D-01 
6  3=  0.1641D-02 
G  5=-0.4545D-06 

G  1=  0,2397D-01 
6  3=  0.1595D-02 
G  5=-0.4731D-06 

6  1=  0.2362D-01 
G  3=  0.1550D-02 
6  5=-0.4919D-06 

G  1=  0,23260-01 
6  3=  0. 15060-02 
G  5=-0. 51080-06 

G  1=  0.2289D-01 
G  3=  0.14620-02 
G  5=-0. 52960-06 

G  1=  0,22510-01 
G  3=  0.14180-02 
G  5=-0. 54840-06 

G  1=  0.22120-01 
G  3=  0.13760-02 
G  5=-0. 56720-06 

G  1=  0,21720-01 
G  3=  0,13330-02 
6  5=-0, 58570-06 

G  1=  0. 21310-01 
G  3=  0.12910-02 
G  5=-0. 60400-06 

G  1=  0.2089D-01 
G  3=  0,12500-02 
G  5=-0,6221D-06 


GP  l=-0, 29480-01 
GP  3=-0. 47060-02 
GP  5=-0, 18000-05 

GP  l=-0, 30630-01 
GP  3=-0. 46600-02 
GP  5=-0, 18240-05 

GP  l=-0, 31780-01 
GP  3=-0. 46140-02 
GP  5=-0. 18440-05 

GP  l=-0. 32940-01 
GP  3=-0, 45670-02 
GP  5=-0. 18610-05 

GP  l=-0. 34080-01 
GP  3=-0. 45190-02 
GP  5=-0. 18740-05 

GP  l=-0. 35220-01 
GP  3=-0, 44710-02 
GP  5=-0, 18820-05 

GP  l=-0, 36350-01 
GP  3=-0, 44220-02 
GP  5=-0, 18850-05 

GP  l=-0. 37470-01 
GP  3=-0, 43730-02 
GP  5=-0, 18840-05 

GP  1 =-0.38570-01 
GP  3=-0. 43230-02 
GP  5=-0, 18760-05 

GP  l=-0, 39650-01 
GP  3=-0. 42730-02 
GP  5=-0, 18640-05 

GP  l=-0. 40700-01 
GP  3=-0, 42240-02 
GP  5=-0. 18440-05 

GP  l=-0, 41740-01 
GP  3=-0, 41730-02 
GP  5=-0, 18190-05 

GP  l=-0. 42740-01 
GP  3=-0, 41230-02 
GP  5=-0. 17860-05 


GPP  l=-0,1152D+00 
GPP  3=  0,45220-02 
GPP  5=-0. 25300-05 

GPP  l=-0. 11540+00 
GPP  3=  0,46140-02 
GPP  5=-0. 22120-05 

GPP  l=-0. 11540+00 
GPP  3=  0.46910-02 
GPP  5=-0, 18590-05 

GPP  l=-0, 11500+00 
GPP  3=  0,47560-02 
GPP  5=-0. 14720-05 

GPP  l=-0. 11440+00 
GPP  3=  0.48110-02 
GPP  5=-0. 10500-05 

GPP  l=-0. 11350+00 
GPP  3=  0.48570-02 
GPP  5=-0. 59280-06 

GPP  l=-0. 11230+00 
GPP  3=  0.48950-02 
GPP  5=-0, 99550-07 

GPP  l=-0. 11080+00 
GPP  3=  0,49290-02 
GPP  5=  0,42930-06 

GPP  l=-0, 10900+00 
GPP  3=  0,49580-02 
GPP  5=  0.99370-06 

GPP  l=-0, 10690+00 
GPP  3=  0,49840-02 
GPP  5=  0,15930-05 

GPP  l=-0. 10450+00 
GPP  3=  0,50090-02 
GPP  5=  0.22270-05 

GPP  l=-0. 10190+00 
GPP  3=  0.50330-02 
GPP  5=  0.28950-05 

GPP  l=-0. 98920-01 
GPP  3=  0,50570-02 
GPP  5=  0.35950-05 


GPPP  l=-0. 391 60-01 
GPPP  3=  0.10070-01 
GPPP  5=  0.30180-04 

GPPP  l=-0. 10300-01 
GPPP  3=  0.84780-02 
GPPP  5=  0.33530-04 

GPPP  1=  0.18800-01 
GPPP  3=  0.71150-02 
GPPP  5=  0.36960-04 

GPPP  1=  0.48040-01 
GPPP  3=  0.59660-02 
GPPP  5=  0.40450-04 

GPPP  1=  0.77360-01 
GPPP  3=  0.50100-02 
GPPP  5=  0.43980-04 

GPPP  1=  0.10670+00 
GPPP  3=  0.42300-02 
GPPP  5=  0.47540-04 

GPPP  1=  0.13600+00 
GPPP  3=  0.36090-02 
GPPP  5=  0,51110-04 

GPPP  1=  0,16510+00 
GPPP  3=  0,31320-02 
GPPP  5=  0.54660-04 

GPPP  1=  0.19410+00 
GPPP  3=  0.27850-02 
GPPP  5=  0.58190-04 

GPPP  1=  0.22300+00 
GPPP  3=  0.25540-02 
GPPP  5=  0.61670-04 

GPPP  1=  0. 25160+00 
GPPP  3=  0.24300-02 
GPPP  5=  0.65070-04 

GPPP  1=  0.27990+00 
GPPP  3=  0.24010-02 
GPPP  5=  0.68380-04 

GPPP  1=  0.30800+00 
GPPP  3=  0.24580-02 
GPPP  5=  0.71570-04 
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X=0,50 


X=0,51 


X=0,53 


X=0,54 


X=0.55 


X=0.56 


X=0.57 


X=0,58 


X=0.59 


X=0.60 


X=0.61 


X=0,62 


G  1=  0.2045D-01  GP  1=-0.4371D-01  GPP  1=-0.9571D-01  GPPP  1=  0.3357D+00 

6  3=  0.1209D-02  GP  3=-0. 407211-02  GPP  3=  0.5082D-02  GPPP  3=  0.25930-02 

G  5=-0.6398D-06  GP  5=-0.1747D-05  GPP  5=  0.4326D-05  GPPP  5=  0.74&0D-04 

G  1=  0,2001D-01  GP  1=-0.4465D-01  GPP  1=-0.9221D-01  GPPP  1=  0.3632D+00 

G  3=  0,11680-02  GP  3=-0.4021D-02  GPP  3=  0.5109D-02  GPPP  3=  0.2797D-02 

G  5=-0.6570D-06  6P  5=-0.1700D-05  GPP  5=  0.5087D-05  GPPP  5=  0.7746D-04 

G  1=  0.1956D-01  GP  1=-0.4556D-01  GPP  1=-0,8844D-01  GPPP  1=  0.3904D+00 

6  3=  0.1128D-02  GP  3=-0,3970D-02  GPP  3=  0,51380-02  GPPP  3=  0,30630-02 

G  5=-0. 67380-06  GP  5=-0.1645D-05  GPP  5=  0,58750-05  GPPP  5=  0.80100-04 

G  1=  0. 19100-01  GP  l=-0. 46420-01  GPP  l=-0. 84400-01  GPPP  1=  0,41720+00 

G  3=  0,10890-02  GP  3=-0.39180-02  GPP  3=  0.51700-02  GPPP  3=  0. 33860-02 

G  5=-0.6899D-06  GP  5=-0, 15820-05  GPP  5=  0,66890-05  GPPP  5=  0.82510-04 

G  1=  0.18630-01  GP  l=-0. 47240-01  GPP  l=-0, 80100-01  GPPP  1=  0,44370+00 

6  3=  0.10500-02  GP  3=-0. 38670-02  GPP  3=  0.52060-02  GPPP  3=  0.37600-02 

G  5=-0, 70540-06  GP  5=-0. 15110-05  GPP  5=  0.75250-05  GPPP  5=  0.84650-04 

G  1=  0,18150-01  GP  l=-0. 48020-01  GPP  l=-0, 75530-01  GPPP  1=  0,46980+00 

G  3=  0.10110-02  GP  3=-0. 38140-02  GPP  3=  0.52450-02  GPPP  3=  0. 41800-02 

G  5=-0.7201D-06  GP  5=-0. 14310-05  GPP  5=  0.83820-05  GPPP  5=  0,86490-04 

G  1=  0,17670-01  GP  l=-0,48750-01  GPP  l=-0, 70700-01  GPPP  1=  0,49570+00 

G  3=  0.97350-03  GP  3=-0. 37620-02  GPP  3=  0.52890-02  GPPP  3=  0,46400-02 

G  5=-0, 73400-06       GP  5=-0. 13430-05  GPP  5=  0,92550-05  GPPP  5=  0.87980-04 

G  1=  0.17180-01  GP  l=-0. 49430-01  GPP  l=-0. 65620-01  GPPP  1=  0,52110+00 

G  3=  0.93610-03       GP  3=-0. 37080-02  GPP  3=  0.53380-02  GPPP  3=  0.51380-02 

G  5=-0. 74700-06  GP  5=-0. 12460-05  GPP  5=  0.10140-04  GPPP  5=  0.89090-04 

G  1=  0.16680-01  GP  l=-0. 50060-01  GPP  l=-0. 60280-01  GPPP  1=  0.54630+00 

G  3=  0,89930-03       GP  3=-0.36550-02  GPP  3=  0.53920-02  GPPP  3=  0. 56680-02 

6  5=-0. 75900-06       GP  5=-0, 11400-05  GPP  5=  0,11040-04  GPPP  5=  0,39780-04 

G  1=  0.16180-01       GP  l=-0, 50640-01  GPP  l=-0, 54690-01  GPPP  1=  0,57110+00 

6  3=  0.86300-03       GP  3=-0. 36010-02  GPP  3=  0.54520-02  GPPP  3=  0.62290-02 

G  5=-0. 76980-06       GP  5=-0. 10260-05  GPP  5=  0.11940-04  GPPP  5=  0.90020-04 

G  1=  0.15670-01       GP  l=-0. 51150-01  GPP  l=-0. 48860-01  GPPP  1=  0.59550+00 

G  3=  0,82730-03       GP  3=-0, 35460-02  GPP  3=  0.55170-02  GPPP  3=  0.68160-02 

G  5=-0. 77950-06       GP  5=-0. 90170-06  GPP  5=  0.12840-04  GPPP  5=  0.89750-04 

G  1=  0.15160-01       GP  l=-0, 51610-01  GPP  l=-0, 42780-01  GPPP  1=  0.61970+00 

G  3=  0.79210-03       GP  3=-0, 34900-02  GPP  3=  0.55880-02  GPPP  3=  0.74280-02 

G  5=-0. 78790-06       GP  5=-0. 76880-06  GPP  5=  0,13730-04  GPPP  5=  0.88930-04 

G  1=  0.14640-01       GP  l=-0. 52010-01  GPP  l=-0. 36460-01  GPPP  1=  0.64350+00 

G  3=  0.75750-03       GP  3=-0, 34340-02  GPP  3=  0,56650-02  GPPP  3=  0. 80620-02 

G  5=-0,7949D-06       GP  5=-0, 62710-06  GPP  5=  0.14620-04  GPPP  5=  0.87520-04 


113 


X=0.63 


X=0.65 


X=0.66 


X=0.67 


X=0.68 


X=0.69 


X=0.70 


X=0,71 


X=0,72 


X=0.73 


X=0.74 


X=0.75 


G  1=  0.1412D-01 
G  3=  0.7235D-03 
G  5=-0.8004D-06 

G  1=  0. 135911-01 
G  3=  0.690011-03 
G  5=-0. 804411-06 

G  1=  0. 130611-01 
G  3=  0,6571D-03 
G  5=-0. 806711-06 

G  1=  0. 125311-01 
G  3=  0.624811-03 
G  5=-0.B074D-06 

G  1=  0. 120011-01 
G  3=  0.5931D-03 
G  5=-0.8062D-06 

G  1=  0.1147D-01 
G  3=  0.562011-03 
G  5=-0.8032D-06 

6  1=  0.109511-01 
G  3=  0.5315D-03 
G  5=-0.7983D-06 

G  1=  0.1042D-01 
G  3=  0.5017D-03 
G  5=-0,7914D-06 

G  1=  0.9892D-02 
G  3=  0.4725D-03 
G  5=-0. 782411-06 

G  1=  0.9370D-02 
G  3=  0.4440D-03 
G  5=-0.7713Ii-06 

G  1=  0.885111-02 
G  3=  0.4162D-03 
G  5=-0.7581D-06 

G  1=  0.833711-02 
G  3=  0.3891D-03 
G  5=-0. 742811-06 

G  1=  0.782911-02 
G  3=  0.3627D-03 
6  5=-0.7253D-06 


6P  l=-0. 523411-01 
GP  3=-0. 337711-02 
GP  5=-0. 476611-06 

GP  1=-0.5261D-01 
GP  3=-0.3319D-02 
GP  5=-0.3175D-06 

GP  1=-0.5280D-01 
GP  3=-0.3260[i-02 
GP  5=-0. 150211-06 

GP  l=-0.5293[i-01 
GP  3=-0.3200D-02 
GP  5=  0. 250111-07 

GP  1=-0.5298D-01 
GP  3=-0.3139D-02 
GP  5=  0.207711-06 

GP  1=-0.5295D-01 
GP  3=-0.3077D-02 
GP  5=  0.3975D-06 

GP  1=-0.5285D-01 
GP  3=-0.3014D-02 
GP  5=  0.5935D-06 

GP  l=-0.5267It-01 
GP  3=-0.2949D-02 
GP  5=  0. 795311-06 

GP  1  =-0.524011-01 
GP  3=-0.2883D-02 
GP  5=  0.1002D-05 

GP  1=-0.5205D-01 
GP  3=-0.2816D-02 
GP  5=  0.1213D-05 

GP  1=-0.5162D-01 
GP  3=-0.2747D-02 
GP  5=  0.1426D-05 

GP  1=-0.5109D-01 
GP  3=-0.2676D-02 
GP  5=  0.1641D-05 

GP  1=-0.504BD-01 
GP  3=-0,2604D-02 
GP  5=  0.1857D-05 


GPP  1=-0.2991D-01 
GPP  3=  0.5749D-02 
GPP  5=  0.1548D-04 

GPP  1=-0.2312D-01 
GPP  3=  0.5B40D-02 
GPP  5=  0.1633D-04 

GPP  1=-0.1611D-01 
GPP  3=  0.5937D-02 
GPP  5=  0.1714D-04 

GPP  l=-0.8863D-02 
GPP  3=  0.6041D-02 
GPP  5=  0.1791D-04 

GPP  l=-0.1394D-02 
GPP  3=  0.615311-02 
GPP  5=  0.1864D-04 

GPP  1=  0.6297D-02 
GPP  3=  0.6271D-02 
GPP  5=  0.193111-04 

GPP  1=  0.1421D-01 
GPP  3=  0.6397D-02 
GPP  5=  0.1991D-04 

GPP  1=  0.2233D-01 
GPP  3=  0.6531D-02 
GPP  5=  0.2044D-04 

GPP  1=  0.3067D-01 
GPP  3=  0.6672D-02 
GPP  5=  0.2089D-04 

GPP  1=  0.3922D-01 
GPP  3=  0.6B21D-02 
GPP  5=  0.2123D-04 

GPP  1=  0,4798D-01 
GPP  3=  0.6977D-02 
GPP  5=  0.2147D-04 

GPP  1=  0. 569511-01 
GPP  3=  0.7142D-02 
GPP  5=  0.2159D-04 

GPP  1=  0.6612D-01 
GPP  3=  0.7314D-02 
GPP  5=  0.215711-04 


GPPP  1=  0.667011+00 
GPPP  3=  0,8716D-02 
GPPP  5=  0.8547D-04 

GPPP  1=  0.6902D+00 
GPPP  3=  0.9389D-02 
GPPP  5=  0,8273D-04 

GPPP  1=  0.7131D+00 
GPPP  3=  0.1008D-01 
GPPP  5=  0.792511-04 

GPPP  1=  0.7357D+00 
GPPP  3=  0.1078D-01 
GPPP  5=  0.7496D-04 

GPPP  1=  0.7580D+00 
GPPP  3=  0.1150D-01 
GPPP  5=  0.6982D-04 

GPPP  1=  0.7800MO0 
GPPP  3=  0.1223D-01 
GPPP  5=  0,6378D-04 

GPPP  1=  0.8018D+00 
GPPP  3=  0.1297D-01 
GPPP  5=  0.567511-04 

GPPP  1=  0.8233D+00 
GPPP  3=  0.1373D-01 
GPPP  5=  0.4870D-04 

GPPP  1=  0.8445D+00 
GPPP  3=  0.1449D-01 
GPPP  5=  0.3955D-04 

GPPP  1=  0.8655D+00 
GPPP  3=  0.1527D-01 
GPPP  5=  0.2923I>-04 

GPPP  1=  0.8863D+00 
GPPP  3=  0.1605D-01 
GPPP  5=  0.1768D-04 

GPPP  1=  0.9068M00 
6PPP  3=  0.1685D-01 
6PPP  5=  0.4B24I)-05 

GPPP  1=  0.9272D+00 
GPPP  3=  0.1765D-01 
GPPP  5=-0.9407D-05 
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X=0.76 


X=0,77 


X=0.78 


X=0.79 


X=0.80 


X=0.81 


X=0.82 


X=0.83 


X=0.84 


X=0,85 


X=0.86 


X=0.87 


X=0.88 


G  1=  0,73280-02 
G  3=  0,3370D-03 
G  5=-0.7057D-06 

6  1=  0.6834D-02 
G  3=  0.3121D-03 
G  5=-0.6839D-06 

G  1=  0,6348D-02 
G  3=  0.2880D-03 
G  5=-0.6600D-06 

G  1=  0,5873D-02 
G  3=  0,2646D-03 
G  5=-0,6340D-06 

G  1=  0.5407D-02 
G  3=  0,24200-03 
G  5=-0, 60610-06 

G  1=  0.49530-02 
G  3=  0,22030-03 
G  5=-0, 57630-06 

G  1=  0.45120-02 
G  3=  0.19940-03 
G  5=-0. 54470-06 

G  1=  0.40840-02 
G  3=  0.17940-03 
G  5=-0. 51140-06 

G  1=  0.36710-02 
G  3=  0.16030-03 
G  5=-0. 47670-06 

G  1=  0.32730-02 
G  3=  0.14220-03 
6  5=-0. 44070-06 

G  1=  0.28930-02 
G  3=  0,12490-03 
G  5=-0.4037D-06 

G  1=  0,25300-02 
G  3=  0.10870-03 
G  5=-0.3659D-06 

6  1=  0.21870-02 
G  3=  0.93430-04 
G  5=-0. 32760-06 


GP  l=-0. 49770-01 
GP  3=-0. 25300-02 
GP  5=  0.2072D-05 

GP  l=-0. 48970-01 
GP  3=-0. 24540-02 
GP  5=  0,22840-05 

GP  l=-0, 48070-01 
GP  3=-0, 23760-02 
GP  5=  0.24920-05 

GP  1 =-0.47070-01 
GP  3=-0. 22960-02 
GP  5=  0. 26940-05 

GP  l=-0. 45970-01 
GP  3=-0, 22140-02 
GP  5=  0.28880-05 

GP  l=-0. 44770-01 
GP  3=-0, 21300-02 
GP  5=  0,30720-05 

GP  l=-0, 43460-01 
GP  3=-0, 20440-02 
GP  5=  0,32420-05 

GP  l=-0. 42050-01 
GP  3=-0. 19550-02 
GP  5=  0.33970-05 

GP  l=-0. 40530-01 
GP  3=-0. 18640-02 
GP  5=  0.35340-05 

6P  l=-0. 38900-01 
GP  3=-0. 17700-02 
GP  5=  0.36500-05 

GP  l=-0. 37160-01 
GP  3=-0. 16740-02 
GP  5=  0.37410-05 

GP  l=-0. 35300-01 
GP  3=-0. 15750-02 
GP  5=  0.38040-05 

GP  l=-0. 33330-01 
GP  3=-0. 14730-02 
GP  5=  0.38360-05 


GPP  1=  0.75490-01 
GPP  3=  0.74950-02 
GPP  5=  0.21400-04 

GPP  1=  0.85060-01 
GPP  3=  0,76830-02 
GPP  5=  0,21070-04 

GPP  1=  0,94840-01 
GPP  3=  0,78800-02 
GPP  5=  0. 20550-04 

GPP  1=  0,10480+00 
GPP  3=  0,80860-02 
GPP  5=  0,19840-04 

GPP  1=  0,11500+00 
GPP  3=  0,82990-02 
GPP  5=  0,18920-04 

GPP  1=  0,12530+00 
GPP  3=  0,85210-02 
GPP  5=  0.17760-04 

GPP  1=  0,13590+00 
GPP  3=  0.87520-02 
GPP  5=  0.16360-04 

GPP  1=  0.14660+00 
GPP  3=  0.89910-02 
GPP  5=  0.14680-04 

GPP  1=  0,15750+00 
GPP  3=  0.92390-02 
GPP  5=  0.12710-04 

GPP  1=  0.16860+00 
GPP  3=  0.94960-02 
GPP  5=  0.10430-04 

GPP  1=  0.17990+00 
GPP  3=  0.97620-02 
GPP  5=  0.78080-05 

GPP  1=  0,19140+00 
GPP  3=  0,10040-01 
GPP  5=  0.48280-05 

GPP  1=  0.20300+00 
GPP  3=  0.10320-01 
GPP  5=  0.14620-05 


GPPP  1=  0.94730+00 
GPPP  3=  0,18460-01 
GPPP  5=-0, 25090-04 

GPPP  1=  0,96720+00 
GPPP  3=  0,19280-01 
GPPP  5=-0, 42300-04 

GPPP  1=  0,98680+00 
GPPP  3=  0,20110-01 
GPPP  5=-0. 61 130-04 

GPPP  1=  0.10060+01 
GPPP  3=  0.20940-01 
GPPP  5=-0. 81640-04 

GPPP  1=  0.10260+01 
GPPP  3=  0,21790-01 
GPPP  5=-0, 10390-03 

GPPP  1=  0,10450+01 
GPPP  3=  0.22640-01 
GPPP  5=-0. 12810-03 

GPPP  1=  0.10640+01 
GPPP  3=  0.23500-01 
GPPP  5=-0. 15420-03 

GPPP  1=  0.10830+01 
GPPP  3=  0,24370-01 
GPPP  5=-0, 18230-03 

GPPP  1=  0.11010+01 
GPPP  3=  0.25240-01 
GPPP  5=-0. 21260-03 

GPPP  1=  0.11200+01 
GPPP  3=  0.26130-01 
GPPP  5=-0. 24510-03 

GPPP  1=  0,11380+01 
GPPP  3=  0. 27020-01 
GPPP  5=-0. 28000-03 

GPPP  1=  0.11560+01 
GPPP  3=  0.27920-01 
GPPP  5=-0. 31730-03 

GPPP  1=  0. 11740+01 
GPPP  3=  0.28830-01 
GPPP  5=-0. 357 10-03 
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X=0.89 


X=0.90 


X=0.91 


X=0.92 


X=0.93 


X=0,94 


X=0.95 


X=0.96 


X=0.97 


X=0,98 


X=0.99 


X=1.00 


G  1=  0.1863D-02 
G  3=  0.7921D-04 
G  5=-0.2892D-06 

G  1=  0.1562D-02 
G  3=  0.6606D-04 
G  5=-0. 251011-06 

G  1=  0.1283D-02 
G  3=  0.5400D-04 
G  5=-0, 213411-06 

G  1=  0.102811-02 
G  3=  0,430611-04 
G  5=-0.1770D-06 

G  1=  0.7979D-03 
G  3=  0.3327D-04 
G  5=-0.1422D-06 

G  1=  0.594411-03 
G  3=  0.2467D-04 
6  5=-0. 109711-06 

G  1=  0.4185D-03 
G  3=  0.1730D-04 
G  5=-0.7992D-07 

6  1=  0.2716D-03 
G  3=  0.1118D-04 
B  5=-0.5370D-07 

G  1=  0.1550D-03 
G  3=  0.6350D-05 
G  5=-0,3175D-07 

G  1=  0.6994D-04 
G  3=  0.2853D-05 
G  5=-0.1488D-07 

G  1=  0.1783D-04 
G  3=  0. 723211-06 
G  5=-0. 397911-08 

G  1=  0.0 
G  3=  0,0 
G  5=  0.0 


GP  1=-0.3124D-01 
GP  3=-0. 136911-02 
GP  5=  0.3832D-05 

GP  l=-0. 290311-01 
GP  3=-0. 126111-02 
GP  5=  0. 378711-05 

GP  1=-0.2670D-01 
GP  3=-0.1150D-02 
GP  5=  0.369911-05 

GP  1=-0.2425D-01 
GP  3=-0.1036D-02 
GP  5=  0.356111-05 

GP  1=-0.2167P-01 
GP  3=-0.9192D-03 
GP  5=  0.3368D-05 

GP  1=-0.1897D-01 
GP  3=-0.7987D-03 
GP  5=  0.3116D-05 

GP  1=-0.1614P-01 
6P  3=-0.6748D-03 
GP  5=  0.2798D-05 

GP  1=-0.1318D-01 
GP  3=-0.5473[i-03 
GP  5=  0.2409D-05 

GP  1=-0.1008D-01 
GP  3=-0,4162D-03 
GP  5=  0.1941D-05 

GP  l=-0. 685811-02 
GP  3=-0.2813D-03 
GP  5=  0.1389D-05 

GP  l=-0.3497D-02 
GP  3=-0.1426D-03 
GP  5=  0.7441D-06 

GP  1=  0.0 
GP  3=  0.0 
GP  5=  0.0 


GPP  1=  0.2149D+00 
GPP  3=  0.1061D-01 
GPP  5=-0.2315D-05 

GPP  1=  0.2269D+00 
GPP  3=  0.1092D-01 
GPP  5=-0. 653111-05 

GPP  1=  0.2391M00 
GPP  3=  0.1123D-01 
GPP  5=-0. 112111-04 

GPP  1=  0.2514D+00 
GPP  3=  0.1155D-01 
GPP  5=-0.1639D-04 

GPP  1=  0.2640D+00 
GPP  3=  0.1188D-01 
GPP  5=-0.2210D-04 

GPP  1=  0.2767D+00 
GPP  3=  0.1222D-01 
GPP  5=-0.2836D-04 

GPP  1=  0.2896D+00 
GPP  3=  0.1257D-01 
GPP  5=-0.3521D-04 

GPP  1=  0.3026D+00 
GPP  3=  0.1293D-01 
GPP  5=-0.4269D-04 

GPP  1=  0.3159D+00 
GPP  3=  0.1330D-01 
GPP  5=-0.5083D-04 

GPP  1=  0,3293D+00 
GPP  3=  0.13680-01 
GPP  5=-0.5966D-04 

GPP  1=  0. 342811+00 
GPP  3=  0.1406D-01 
GPP  5=-0, 692311-04 

GPP  1=  0.3566D+00 
GPP  3=  0,1446D-01 
GPP  5=-0. 795811-04 


GPPP  1=  0.1192D+01 
GPPP  3=  0.297511-01 
GPPP  5=-0.3996D-03 

GPPP  1=  0. 121011+01 
GPPP  3=  0.3O67D-O1 
GPPP  5=-0.4449D-03 

GPPP  1=  0.1228M01 
GPPP  3=  0.3161D-01 
GPPP  5=-0.4930D-03 

GPPP  1=  0.1245D+01 
GPPP  3=  0.3256D-01 
GPPP  5=-0.5441D-03 

GPPP  1=  0.1263D+01 
GPPP  3=  0. 335111-01 
GPPP  5=-0.5983D-03 

GPPP  1=  0.1280D+01 
GPPP  3=  0.3447D-01 
GPPP  5=-0,6558D-03 

GPPP  1=  0.1297D+01 
GPPP  3=  0.3545D-01 
GPPP  5=-0.7165D-03 

GPPP  1=  0.1315D+01 
GPPP  3=  0.3643D-01 
GPPP  5=-0.7808D-03 

GPPP  1=  0.133211+01 
GPPf  3=  0.3743D-01 
GPPP  5=-0. 848711-03 

GPPP  1=  0.1349D+01 
GPPP  3=  0.3843D-01 
GPPP  5=-0. 920311-03 

GPPP  1=  0.1366D+01 
GPPP  3=  0,3945D-01 
GPPP  5=-0. 995811-03 

GPPP  1=  0.1383D+01 
GPPP  3=  0.4045D-01 
GPPP  5=-0.1067D-02 


IER=  0 
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DIHENSIONLESS  STREAh 

X=0.0 


X=0,10 


X=0.20 


X=0.30 


HU=0.0 

PSI= 

0.25000ri-0i 

MU=0.10 

PSI= 

0. 2499711-01 

HU=0.20 

PSI= 

0. 2496011-01 

hU=0,30 

PSI= 

0.2479711-01 

Mil- A  AT. 

rSl- 

0. 2436011-01 

HU=0.50 

PSI= 

0.2343711-01 

MU=0.60 

PSI= 

0,2176011-01 

MU=0.70 

PSI= 

0.1899711-Oi 

MU=0,80 

PSI= 

0. 1476011-01 

HLI=0.90 

PSI= 

0.8597511-02 

HU=1.00 

PSI= 

0.0 

HU=0.0 

PSI= 

0.2499411-01 

«U=0.10 

PSI- 

0.24974[i-Ci 

«U=0.20 

PSI- 

0.24886D-01 

HLl=0.30 

PSI= 

0,2464611-01 

rtLl=0i40 

PSI= 

0.2411611-01 

hLI=0,50 

PSI= 

0.2309''D-C1 

HU=0.60 

PSI= 

0.2134611-01 

HU=0.70 

PSI= 

0.18550D-01 

HU=0,80 

PSI= 

0.143481.-01 

hU=0.50 

PSI= 

0.8322411-02 

HU=1.00 

PSI= 

0,0 

hU=0,0 

PSI= 

0.24767D-01 

HU=0.10 

PSI= 

0.24715ri-01 

MU=0.20 

PSI= 

0.2453311-01 

HU=0,30 

PSI= 

0.2415211-01 

MU=0.40 

PSI= 

0,2345311-01 

«U=0.50 

PS1= 

0.22270D-01 

KU=0,60 

PSI= 

0.20391P-01 

HU=0.70 

PSI= 

0. 1755511-01 

HU=0.80 

P3I= 

0. 1345611-01 

HU=0.90 

PSI= 

0. 7738011-02 

HU=hOO 

PSI= 

0.0 

hU=0,0 

PSI= 

0. 2384411-01 

HU=0.10 

PSI  = 

0. 2376311-01 

HU=0,20 

PSI= 

0.23501D-01 

HU=0.30 

PSI= 

0. 2300111-01 

MU=0.40 

PSI= 

0.2216711-01 

HLI=0.50 

PSI= 

0.2086511-01 

MU=0.60 

PSI= 

0.1892411-01 

«U=0,70 

PSI= 

0.1613411-01 

HU=0.80 

PSI= 

0.1224711-01 

HU=0,90 

PSI= 

0. 6976711-02 

HU=1.00 

PSI= 

0.0 

ION  PSKXjMU) 
X=0.40 


X=0,50 


X=0,60 


.X=0,70 


Hll=0.0 

PSI= 

0.2184011-01 

HU=0.10 

PSI= 

0.21743I.-01 

HU=0.20 

PSI= 

0.2143911-01 

MU=0.30 

PSI= 

0.20883D-01 

0.7000011-01 

MU=0.50 

PSI= 

0.166891i-0i 

«U=0,60 

PSI= 

0. 1681411-01 

HU=0.70 

PSI= 

0.14214D-01 

HU=0.80 

PSI= 

0.1069811-Oi 

MLI=0.90 

PSI= 

0. 6043011-02 

HU=1.00 

P3I= 

0.0 

HU=0,0 

PSI= 

0.1863911-01 

hll=0.10 

PSI= 

0,1 854211-0  i 

HU=0.20 

PSI= 

0.1824211-01 

HLi=0.30 

PSI= 

0.1770511-01 

HU-O.tO 

rSi  = 

0,1687711-01 

HU=0.50 

P3I= 

0.1568111-Oi 

HU=0.60 

PSI= 

0,14019[t-01 

MU=0.70 

PSI  = 

0.1i772D-01 

MLl=0,80 

PSI= 

0.8798511-02 

MU=0.90 

PSI= 

0.4935811-02 

MU=1,00 

P3I= 

0.0 

MU=0.0 

PSI= 

0.14427D-01 

HU=0.10 

PSI= 

0,14345D-01 

HU=0,20 

PSI= 

0. 1408911-01 

HLi=0.30 

PSI= 

0,13638iK'i 

MLI=0.40 

PSI= 

0.12955D-01 

HU=0.50 

PSI= 

C.  11 98611-0 1 

HU=0.60 

PSI= 

0.1066511-01 

MLl=0.70 

PSI= 

0,8909811-02 

«U=0,80 

PSI= 

0,6623611-02 

MU=0.90 

PSI= 

0.3695411-02 

HU=1.00 

P3I= 

0.0 

HU=0,0 

PSI= 

0.9663811-02 

HU=0.10 

PSI= 

0,96046D-02 

«U=0.20 

PSI= 

0,9422511-02 

Mll=0.30 

PSI= 

0,9103711-0"' 

«U=0.40 

PSI= 

0. 8625411-02 

«U=0.50 

PSI= 

0.7955811-02 

hLI=0.60 

PSI= 

0.7054011-02 

MLI=0.70 

PSI  = 

0.5870211-02 

MU=0,80 

PSI= 

0. 4346111-02 

MU=0.90 

PSI= 

0.2414511-02 

HU=1.00 

PSI= 

0,0 

X=0,30 


X^l.OO 


X=0.90 


HU=0.0 

PSI= 

0.30429D-02 

HU=0,0 

PSI= 

0.0 

(iLl=0.10 

P3I= 

0.5010611-02 

«U=0.10 

P3I= 

0,0 

MU=0.20 

PSI= 

0.49114D-0'' 

MU=0.20 

V  ♦  V 

HU=0.30 

PSI  = 

0,4733811-02 

HU=0,30 

PSI= 

0.0 

«U=0,40 

PSI= 

0.44316D-02 

HU=0.40 

PSI= 

0.0 

«U=0.50 

P3I= 

0,41244[i-O2 

HU=0,50 

PSI= 

0.0 

i1U=0,60 

PSI= 

0,36474D-02 

HU=0.60 

PSI= 

0,0 

HU=0.70 

PSI= 

0.30266D-02 

riU=0,70 

PSI= 

0.0 

HU=0,80 

PSI= 

0.22338D-02 

MU=0.80 

P3I= 

0,0 

HLI=0.90 

PSI= 

0.12370D-02 

MU=0,?0 

PSI= 

0.0 

«U=1.00 

PSI= 

0,0 

HU=1.00 

PSI:: 

0,0 

HU=O,0 

PSI= 

0.14423D-02 

)1U=0.10 

PSI= 

0.1452611-02 

HU=0,20 

PSI= 

0.14230D-02 

HIJ=0,30 

PSI  = 

0.1371 711-02 

HU=0.40 

PSI= 

0.12956D-02 

MU=0,50 

PSI= 

0,1190411-02 

HU=0,60 

PSI= 

O.i  050711-0  2 

MU=0.70 

PSi  = 

0 , 8700011-03 

HU=0,80 

PSI= 

0.6406311-03 

I1U=0,90 

PSI  = 

0,3538611-03 

HU=1,00 

PSI= 

0.0 
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X=0.10 


i'lblKitlui  iUN  ViUflL..' \miUj 

X=0,40 

HU=0.0 

y= 

0.0 

MU=0.0 

0.65226D-0i 

MU=0.10 

W= 

0. 4000011-02 

hU=0.10 

0.65680ri-Ci 

HU=0,20 

H= 

0.1 600011-0  i 

rtU=0.20 

w= 

0.6716611-01 

HLI=0,30 

W= 

0,360C0[i-01 

hU=y.30 

y= 

0.70042D-01 

hLI=0.40 

II. 

w= 

y.  6400011-01 

HU=0.40 

w= 

A     TAATCT.  A-l 

0,7483511-01 

nu-D  t  M 

u- 
W- 

A   1  AAAAriiAA 

y .  lyyuyiityu 

Ul  l*.A  CTA 

0,c2118D-yl 

MLI=0,60 

H= 

0.1440011+00 

hU=0.60 

W= 

0.9236911-01 

H\\=().lCi 

ft" 

0.1960011+00 

Mil-A  7A 

nu-v./y 

U- 

y.lODbWU+OO 

I lU    V  ♦  UV 

u= 

0. ''540011+ 00 

UtI-A  OA 

u- 
tt- 

y.li:2//ii+00 

nu-y  t TV 

u- 

hU=0.90 

W= 

0. 1430111+00 

Ml  1=1  .fir; 

ij= 
w- 

V  t  tvvvUi'Tv  -,' 

X=0,50 

hU=1.00 

H= 

0.1665311+00 

UII— A  A 

nu-OlO 

Ll  — 

MU=0,0 

W= 

0.7521511-01 

MU=0.10 

w= 

0.10431Ii-yl 

HU=0,10 

W= 

0.752701.-01 

rtU=0 . 20 

w= 

0.2j239Ii-01 

hU=0,20 

W= 

0.7547611-01 

HU=0.30 

w= 

0,42336Ii-0l 

MU=0.30 

W= 

0. 7595411-01 

MU=0.40 

w= 

0,6811111-01 

HU=0.40 

w= 

0. 7691211-01 

t1Ll=0.jO 

H= 

0. 1008111+00 

KU=0.50 

w= 

0.7864211-01 

HU=0.60 

W= 

0.1 405611+00 

HU=0.60 

w= 

0.6149811-01 

nU=0 1 70 

1(1= 

0.1S74JD+00 

HL1=0,70 

w= 

0. 8535311-01 

UM-A  OA 

H- 

0.241411.1+00 

HU=0>80 

w= 

0.9204111-01 

mt-ri  Oi''i 

nu-u I yo 

u- 
H- 

y.i02j6Ii+00 

«U=0.90 

H= 

0.10030I'+00 

HI  1-1  An 

11- 

A  T7AQ7T(1AA 

X=0.60 

HU=1.00 

W= 

0.1107811+00 

u- 

v.ii/iii'  yi 

HU=0.0 

W= 

0.7639711-01 

Hibfi  1A 
n\j— V » 1 V 

u- 

A  'l7='.Aijri-A1 

v.t/jyyi'  yi 

HU=0,10 

0,76272Ei-Cl 

u- 

y.ojotJoii  yi 

HU=0.20 

w= 

0,7590111-01 

Mil- A  1l\ 

Ll- 

W- 

0.4V8ri.'li-01 

HU=0.30 

w= 

0.75287ri-Cl 

u~ 

H- 

0,6956311-01 

HU=0,40 

w= 

0.7444411-01 

nU=0.50 

W= 

0.94912D-01 

ML1=0,50 , 

H= 

0.7339711-01 

hU=0.60 

U= 

0.1259411+00 

«U=0,60 

0.72191D-01 

MU=0.70 

w= 

0. 1626311+00 

I1U=0.70 

0.7089711-01 

MU=0.80 

H= 

0.2050011+00 

H1J=0,80 

0.6961311-01 

tin     A.     r\  A 

(1U=0.90 

W= 

0,2530111+00 

hU=0.90 

w= 

0.68466B-01 

nU-iiOO 

tl_ 

0. 3066811+00 

X=0.70 

HU=1.00 

w= 

0,6761211-01 

Mll-A  A 

u- 

y.i/Jiioii-oi 

HU=0.0 

w= 

0,6892011-01 

Mll-A  1A 

u- 
w- 

y .iDj/ui  01 

Hli=0.10 

w= 

0,6874211-01 

Mll-A  ICi 

nu-y  t  iV 

u- 

A  cio^nn  Ai 
y .  ji-JO/U  01 

hU=0.20 
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